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ABSTRACT 


A  method  for  computing  three-dimensional  flow  over  an  ogival  body  at  angle  of 
attack  is  described.  An  approximate  set  of  governing  equations  is  derived  for 
viscous  flows  which  have  a  primary  flow  direction.  The  derivation  is  done  in  a 
coordinate  independent  manner,  and  the  resulting  equations  are  expressed  in  terms 
of  tensors.  In  keeping  with  the  inherent  generality  of  the  tensor  formulation,  a 
two-level  second-order  accurate  marching  procedure  is  derived  for  general  tensor¬ 
like  equations.  With  this  procedure,  a  three -dim Msional  turbulent  flow  can  be 
solved  in  any  coordinate  system  by  marching  along  the  assumed  primary  flow  direction. 
General  tube-like  coordinates  are  developed  for  a  class  of  geometries  applicable 
to  flows  between  tubular  surfaces.  The  coordinates  are  then  particularized  to  the 
flow  field  bounded  between  an  ogival  body  at  angle  of  attack  and  its  bow  shock. 
Unlike  the  ogival  body  surface,  the  bow  shock  surface  is  not  known  in  advance  of  the 
solution  but  instead  must  be  computed  as  the  solution  develops.  One  marching  step 
of  the  solution  process  is  broken  down  into  several  steps.  First,  the  bow  shock 
surface  is  discretely  extended  by  an  iteration  of  explicit  local  solutions.  The 
bow  shock  surface  is  then  smoothly  extended  to  provide  a  best  fit  to  the  discrete 
shock  data.  Tube-like  coordinates  are  generated  and  finally  the  second  order 
numerical  scheme  is  applied  to  advance  the  fully  viscous  solution  to  the  next 

station.  / 
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A  Method  for  Computing  Three-Dimensional  Viscous  Flows 
Over  an  Ogival  Body  at  Angle  of  Attack 


INTRODUCTION 


An  important  consideration  in  the  design  of  supersonic  missiles  is  the  accurate 
prediction  of  both  the  pressure  distribution  and  heat  transfer  loads  about  the  body. 
Although  the  combination  of  inviscid  flow  theory  and  three-dimensional  boundary 
layer  theory  may  be  adequate  to  predict  the  flows  about  ogival  bodies  at  small 
angles  of  attack,  these  analyses  used  separately  are  inadequate  at  larger  angles 
of  attack.  At  larger  angles  of  attack,  a  strong  viscous-inviscid  interaction  occurs 
on  the  lee  side  of  the  body  leading  to  the  formation  of  a  pair  of  vortices  symmetric 
about  the  lee  body  generator  and  an  accurate  flow  field  prediction  under  these  con¬ 
ditions  requires  a  solution  which  considers  the  mutual  interaction  between  the 
viscous  layer  and  the  nominally  inviscid  flow.  Since  increases  in  the  body  lift 
to  drag  ratio  and  in  local  heat  transfer  rate  on  the  lee  side  of  the  body  are 
associated  with  the  formation  of  the  vortices,  an  accurate  method  of  predicting 
the  lee  side  interacting  flow  field  is  necessary  to  insure  both  the  effective 
operation  and  structural  integrity  of  supersonic  vehicles. 

Successful  predictions  of  the  flow  about  ogival  bodies  at  angle  of  attack 
require  an  accurate  flow  model  in  which  the  strong  viscous-inviscid  coupling  on  the 
lee  side  of  the  vehicle  at  large  angle  of  attack  is  modeled  correctly.  Boundary 
layer  theory  is  not  adequate  to  describe  the  development  of  viscous  flow  phenomena 
as  complex  as  vortex  rollup  since  the  basic  assumption  that  the  boundary  layer  ie 
thin  compared  to  a  typical  body  scale  is  invalid.  Thus,  a  three-dimensional 
interacting  boundary  layer  theory  would  be  inadequate  for  description  of  the  problem 
Although  a  numerical  solution  to  the  full  Navier-Stokes  equations  would  provide  a 
theory  with  tins  generality  tc  aucc&sp  fully  predict  such  complex  flows, 

the  required  computer  time  and  storage  indicate  that  three-dimensional  Navier- 
Stokes  solutions  should  be  used  only  if  no  suitable  alternative  exists.  An  optimum 
analysis  would  possess  the  general  three-dimensional  nature  of  the  Navier-Stokes 
equations  but  would  not  be  limited  by  the  large  running  time  and  storage  require¬ 
ments  associated  with  three-dimensional  Navier-Stokes  solutions. 

The  problem  of  predicting  the  flow  field  about  sharp  nosed  ogival  bodies  at 
incidence  has  been  under  investigation  for  the  designer  in  his  consideration  of  aero 
dynamic  forces  and  heating  loads.  At  supersonic  speeds  the  flow  over  ogival  bodies 
at  incidence  may  be  thought  of  as  having  a  component  aligned  with  the  free -stream 
flow  direction,  which  is  little  affected  by  viscous  forces,  and  cross  flow  com¬ 
ponents  flowing  around  the  body,  which  can,  at  large  angles  of  attack,  become 
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viscous  dominated.  This  picture  of  ogival  flow  fields  is  borne  out  by  various 
experimental  investigations  (see  Refs.  1,  2,  3,  and  4).  These  investigations 
have  found  that  at  zero  angle  of  attack  a  symmetric  flow  field  develops  about  the 
body,  in  which  no  cross  flows  exist.  As  the  angle  of  attack  is  increased  from 
zero,  a  cross  flow  pattern  begins  to  develop  proceeding  from  the  windward  symmetry 
plane  toward  the  leeward  symmetry  plane .  When  the  angle  of  attack  becomes  large 
enough,  the  cross  flow  adverse  pressure  gradient  existing  near  the  lee  symmetry 
plane  becomes  sufficiently  strong  and  the  cross  flow  separates,  leading  to  the 
development  of  a  pair  of  vortices  symmetric  about  the  lee  symmetry  plane  as  depicted 
in  Fig.  1.  These  vortices  may  or  may  not  have  imbedded  shocks  associated  with 
them.  The  flow  as  depicted  in  Fig.  1  has  been  verified  for  both  laminar  and  tur¬ 
bulent  flows  for  speeds  in  the  supersonic  and  hypersonic  regimes.  The  available 
experimental  data  indicates  that  flow  separation  (reversal)  in  the  axial  direction 
usually  is  not  associated  with  the  vortex  development. 

The  majority  of  previous  attempts  to  predict  flows  of  this  type  may  be  cate¬ 
gorized  as  either  solutions  of  the  inviscid  flow  equations  or  solutions  to  the 
three-dimensional  boundary  layer  equations.  Even  with  the  simplifying  assumptions 
of  these  two  approaches  the  equations  must  still  be  solved  numerically.  Further¬ 
more,  most  previous  investigations  have  been  limited  to  the  problem  of  flow  over 
conical  bodies  in  which  axial  invariance  is  assumed  in  order  to  eliminate  deriva¬ 
tives  in  the  direction  of  the  cone  axis.  Such  analyses  cannot  be  used  to  analyze 
general  body  shapes  or  to  accurately  predict  the  development  of  the  lee  side  vortices 

The  most  common  inviscid  procedure  currently  in  use  involves  numerical  solution 
of  the  time-dependent  inviscid  flow  equations.  The  steady  inviscid  flow  solution  is 
approached  asymptotically  for  large  values  of  time.  The  solutions  of  Moretti 
(Refs.  5  and  6)  are  examples  of  this  type  of  approach.  Other  available  inviscid 
techniques  Include  the  inverse  method  (Ref.  7)  and  the  method  of  integral  relation^ 
(Ref.  8).  MacCormack  and  Warning  (Ref.  9)  have  recently  surveyed  the  available 
inviscid  computational  procedures.  An  obvious  disadvantage  of  a  purely  irvieclJ 
solution  to  this  problem  is  the  failure  to  account  for  viscous  effects.  Viscous 
forces  may  be  accounted  for,  however,  by  making  use  of  the  inviscid  pressure  distri¬ 
bution  to  solve  the  three-dimensional  boundary  layer  equations  (see  Ref.  10). 

This  procedure  can  give  accurate  predictions  of  ogival  flow  fields  provided  the 
angle  of  attack  is  small  enough  to  prevent  cross  flow  separation  and  hence  does  not 
permit  lee  side  vortices.  However,  the  small  angle  of  attack  constraint  places  a 
severe  restriction  on  the  boundary  layer  type  procedures. 

Because  of  their  complexity,  and  particularly  the  interaction  which  occurs 
between  primary  and  secondary  flows  and  viscous  and  inviscid  regions,  three-dimen¬ 
sional  flows  over  ogival  bodies  at  nontrivial  angles  of  attack  have  been  extremely 
difficult  to  analyze.  There  are  considerable  difficulties  associated  with  the 
synthesis  of  inviscid  flow  analysis  and  boundary  layer  theory  into  a  cohesive  method 
of  analysis.  Among  these  difficulties  are  the  lack  of  applicability  of 
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three-dimensional  boundary  layer  theory,  a  means  for  patching  or  interfacing  bound¬ 
ary  layer  and  rotational  inviscid  flow  regions,  and  the  treatment  of  interaction 
between  viscous  and  inviscid  flow  regions. 

In  efforts  to  develop  methods  for  dealing  with  problems  of  this  type,  Patankar 
&  Spalding  (Ref.  ll),  Caretto,  Curr,  &  Spalding  (Ref.  12),  and  Briley  (Ref.  13) 
devised  numerical  methods  for  solving  approximate  governing  equations  which  are  a 
more  or  less  natural  generalization  of  three-dimensional  boundary  layer  theory.  In 
these  studies,  solutions  were  computed  for  laminar  incompressible  flow  in  straight 
ducts  with  rectangular  cross  sections.  The  governing  equations  were  solved  by  inte¬ 
grating  in  a  primary  flow  coordinate  direction  while  retaining  viscous  stresses  in 
both  transverse  coordinate  directions  as  opposed  to  only  one  direction  for  three- 
dimensional  boundary  layer  theory.  In  addition,  certain  assumptions  were  made  about 
the  behavior  of  pressure  gradient  terms  for  incompressible  flow  to  permit  solution 
by  forward  marching  integration.  Subsequently,  this  general  approach  has  been  used 
to  compute  laminar  incompressible  flow  in  helical  tubes  by  Patankar,  Pratap,  & 

Spalding  (Ref.  14).  A  predictor/corrector  solution  procedure  has  been  developed  by 
Lin  and  Rubin  (Refs.  15  and  l6)  to  solve  the  parabolized  three  dimensional  compres¬ 
sible  Navier  Stokes  equations.  The  numerical  technique  is  implicit  in  one  transverse 
direction  and  iterative  in  the  other.  Helliwell  and  Lubard  (Ref.  17),  Rakich  and  Lubard 
(Ref.  18) ,  .and  Lin  and  Rubin  (Ref.  16)  have  applied  this  method  to  the  problem  of 
flow  over  both  sharp  and  spherically  blunted  cones  at  angle  of  attack. 

Recently  in  companion  studies,  Briley  &  McDonald  (Ref.  19)  anc  McDonald  &  Briley 
(Ref.  20)  have  developed  stable  and  efficient  noniterative  implicit  numerical  tech¬ 
niques  for  application  to  systems  of  coupled  nonlinear  multidimensional  nonelliptic 
equations.  These  general  techniques  were  applied  by  McDonald  &  Briley  (Ref.  20)  to 
the  computation  .by  forward  marching  integration  of  laminar  supersonic  flow  in  rec¬ 
tangular  jets.  Subsequently,  the  laminar  incompressible  straight-duct  analysis  of 
Briley  (Ref.  13)  and  the  improved  numerical  techniques  of  McDonald  &  Briley  (Ref.  20) 
for  compressible  flows  were  extended  and  synthesized  by  Briley  &  McDonald  (Ref.  21) 
and  Eiseman,  McDonald,  and  Briley  (Ref.  22)  into  a  method  for  computing  subsonic 
turbulent  flow  in  curved  ducts.  The  present  study  represents  a  further  generali¬ 
zation  wf  the  latter  method,  to  encompass  graved  art?  highly  me*  lex 

geometries. 

In  the  ogival  body  problem  the  basic  geometry  is  determined  by  both  the  ogival  body  it¬ 
self  and  by  the  bow  shock  propagated  from  the  tip  of  the  body .  Unlike  the  body  shape,  the  oow 
shock  is  not  known  in  advance,  but  must  determined  as  part  of  the  solution.  Since  the  region 
neares  b  the  shock  is  dominated  by  convective  forces ,  and  the  shock  is  treated  as  a  discontin¬ 
uity,  it  is  quite  sufficient  to  perform  a  local  inviscid  analysis  in  that  region;  and,  thereby, 
to  determine  the  shock  location  one  step  in  advance  of  the  fully  viscous  solution  which  is 
being  marched  along  the  axis  of  the  ogival  body.  The  shock  location  is  calculated  numer¬ 
ically  in  terms  of  local  extensions  of  the  existing  coordinate  system.  In  this  way 
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each  new  point  of  shock  location  can  be  explicitly  computed  independent  of 
neighboring  transverse  points;  and,  hence,  remove  any  need  to  construct  a  poten¬ 
tially  costly  global  coordinate  extension  which  would  utlimate'ly  be  discarded. 
Specifically,  for  each  new  shock  point  we  have  a  disti.  ctly  tailored  coordinate 
system  which  is  an  extension  of  the  previous  coordinate  system  that  encompasses 
only  the  explicit  difference  molecule  for  the  point  in  question.  A  distinct  advan¬ 
tage  here  is  that  the  metric  data  necessary  to  write  the  appropriate  equations 
must  only  be  determined  at  the  shock  point  in  question  since  the  metric  data  at 
the  explicit  level  is  already  in  storage  from  the  previous  step.  After  the  shock 
data  has  been  determined  one  is  left  with  a  loop  of  points  one  station  ahead  of  a 
smooth  shock  surface.  The  problem  now  is  to  extend  the  shock  surface  in  a  suffi¬ 
ciently  smooth  and  uniform  manner  so  that  suitable  coordinates  can  be  generated 
for  the  advancement  of  the  fully  viscous  solution  to  the  same  axial  location  of  the 
just  completed  shock  calculation.  The  loop  of  shock  points  may  lack  a  certain  amount 
of  uniformity  by  having  oscillations  in  curvature  if  a  strict  interpolation  were  to 
be  performed.  The  lack  of  uniformity  can  easily  arise  from  the  numerical  nature  of  the 
calculation.  Thus  a  least- squares  spline  procedure  is  employed  to  remove  the  noisy 
oscillations  of  the  data  and  to  produce  a  loop  which  is  smooth  enough  to  possess 
continuous  third  derivatives  (needed  for  a  viscous  calculation)  and  uniform  enough 
to  have  curvature  which  reflects  the  global  structure  of  the  shock  surface.  At 
this  point  one  has  a  well-defined  loop  in  front  of  a  well-defined  surface.  The  sur¬ 
face  is  now  easily  extended  in  a  patchwise  fashion  by  smoothly  joining  polynominal 
surface  elements.  From  here  the  coordinates  are  generated  by  a  linear  deformation 
of  the  ogival  body  surface  into  the  bow  shock  surface.  Tne  rate  of  le formation  is 
controlled  by  a  choice  of  the  linear  deforming  parameter .  This  parameter  is  con- 
strained  to  monot  jTii^ally  and  smoothly  vary  ffcob.  ee*c  to  unity.  A*  this  vmH  at  ton  occurs,  a 
family  of  coordinate  surfaces  is  generated  from  the  ogival  body  to  the  bow  shock. 

When  the  parameter  is  chosen  to  be  a  function  of  only  the  deforming  direction,  the 
deformation  is  uniform  over  the  body  surface.  This  is  ideal  for  global  resolutions 
over  or  within  an  object.  If,  in  addition,  the  deformation  parameter  should  also 
depend  upon  surface  location,  then  local  resolutions  would  be  possible.  In  the 
ogival  body  problem,  however,  the  global  resolutions  are  sufficient  to  adequately 
resolve  the  large  velocity  gradients  associated  with  the  attached  boundary  layer . 

Throughout  the  discussion  to  follow  the  material  will  be  ordered  from  the  more 
general  to  the  more  specific.  In  so  doing  one  has  more  flexibility  and  generality 
at  hand  to  apply  to  the  ogival  body  problem.  Thus  it  is  best  to  start  off  with  a 
discussion  on  the  rationale  behind  the  use  of  curvilinear  coordinate  systems  for 
fluid  dynamic  problems.  While  the  application  -fvt  ridUB  ecordlubte  sjrstPMt  i*  a 
new  idea,  the  coordinate  independent  concept  associated  with  analyses  based  upon 
the  metric  tensor  is  of  greater  utility.  The  metric  tensor,  as  developed  in  the 
initial  discussion,  is  used  to  compare  coordinate  systems  of  varying  degrees  of 
generality.  In  the  next  section  use  is  made  of  the  metric  tensor  to  obtain  a  tensor 
(coordinate  independent)  form  of  the  Navier-Stokes  Equations  which  are  then  approxi¬ 
mated  to  produce  an  initial  value  problem.  The  approximation  is  obtained  from  a 
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neglect  of  viscous  stresses,  and  henc'j  diffusion,  in  an  assumed  primary  flow  direction. 
The  approximation  can  be  viewed  as  a  generalization  of  bound aiy  layer  theory.  The 
primary  flow  direction  is  assumed  to  be  given  by  same  smooth  vector  field.  Since 
the  specification  of  any  vector  field  is  independent  of  coordinates,  the  approxi¬ 
mation  of  the  tensor  form  of  the  Navier-Stokes  equations  is  also  independent  of 
coordinates.  As  a  matter  of  convenience,  however,  the  primary  vector  field  is 
often  chosen  to  be  the  vector  field  associated  with  a  given  coordinate  direction 
of  a  given  system  of  coordinates.  For  the  ogival  body  problem  this  is  done.  To 
maintain  the  generality  of  the  discussion  one  next  considers  the  general  numericr' 
method  rather  than  the  details  of  coordinate  generation.  This  is  consistent  with 
this  general  methodology  since  the  computer  code  is  constructed  in  a  modular  fashion 
in  terms  of  the  metric  data.  After  this  discussion  one  considers  the  construction 
of  coordinates  that  are  suitable  for  the  ogival  body  problems;  however,  some  gen¬ 
erality  is  still  maintained  by  considering  coordinate  systems  that  are  generated 
from  any  two  concentric  tubes.  The  metric  data  is  then  obtained  so  that  the 
arbitrary  two  tube  problem  is  fully  specified  once  the  tubes  are  specified.  In 
the  present  study  the  inner  tube  is  taken  to  be  the  ogival  body  and  the  outer  tube 
is  taken  to  be  the  bow  shock.  One  finally  considers  the  numerical  generation  of 
bow  shock  data,  and  the  use  of  that  data  to  suitably  extend  the  bow  shock  surface 
so  that  a  well-defined  outer  tube  is  specified. 
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THE  METRIC  TENSOR  AND  ITS  RELATIONSHIP  TO  SYSTEMS  OF  COORDINATES 


The  governing  equations  for  a  viscous  fluid  will  be  expressed  and  approximated  in 
generalized  coordinates.  Like  any  physical  process,  the  dynamics  of  a  fluid  is  inde¬ 
pendent  of  coordinates;  and  is,  therefore,  describable  in  terms  of  arbitrary  coordinates. 
The  practical  implication  of  this  coordinate  independence  is  that  the  analyst  has  the 
freedom  to  select  coordinate  systems  which  are  easy  to  construct  and  which  simplify 
the  solution  process. 

In  the  numerical  solution  of  fluid  dynamic  problems  there  are  many  advantages  to 
be  gained  by  judicious  choice  of  coordinates.  The  most  obvious  advantage  is  that  the 
physical  boundaries  of  a  flow  region  can  be  represented  by  coordinate  surfaces.  This 
removes  the  need  for  fractional  cells  and  hence  removes  the  complications 
and  loss  of  accuracy  associated  with  an  interpolation  algorithms  for  the  boundaries. 
Another  advantage  in  the  use  of  generalized  coordinates  is  that  a  uniform  numerical 
method  can  be  used.  The  solution  can  then  be  performed  with  a  fixed  number  of  cells 
in  any  given  direction  and  with  a  uniform  mesh  spacing.  The  result  is  a  simplification 
of  the  computer  logic;  and  hence,  a  savings  in  time  for  both  the  computer  and  the 
programmer.  For  the  ADI  method  of  solution  that  we  use,  the  one-dimensional  rows 
and  columns  each  have  fixed  lengths;  and  hence,  we  are  not  faced  with  the  combinatorical 
problem  of  monitoring  the  lengths  of  rows  and  columns  which  would  otherwise  be  caused 
by  geometric  changes  in  the  boundaries.  In  addition  to  the  above  there  is  the  advantage 
of  an  implicit  mesh  distribution.  The  uniform  mesh  of  computational  space  is  simply 
mapped  into  an  appropriately  distributed  mesh  in  physical  space.  Thus,  the  coordinate 
transformation  can  be  chosen  to  contain  the  distributional  information  as  well  as  the 
boundary  specifications.  The  resolution  of  a  rapidly  changing  solution  is  the  major 
objective  in  selection  of  a  coordinate  mesh  distribution.  A  classical  example  is  the 
resolution  of  attached  boundary  layers  where  the  solution  is  known  to  have  large, 
velocity  gradients.  Another  more  subtle  example  is  the  resolution  of  large  gradients 
in  computational  coordinates  due  to  regions  of  high  curvature  on  the  bounding  surfaces. 
When  the  transformation  contains  the  distributional  information  there  is  no  need.to 
construct  the  apparatus  for  the  discrete  approximation  of  derivatives  on  a  nonuniform  . 
mesh.  This  is  a  savings  in  both  computer  logic  and  storage.  With  a  different  motivation, 
however,  it  may  be  desired  to  automate  the  difference  molecule  so  that  the  numerical 
technique  can  be  changed  with  a  few  parameters.  Changes,  in  practice,  usually  amount 
to  a  selection  between  forward,  backward  or  central  differences.  For  any  given  direction 
we  need  three  parameters  for  first  derivatives  and  three  parameters  for  second  derivatives. 
Thus,  even  with  such  an  automation  of  the  numerical  method  we  save  on  computer  logic 
and  storage.  Specifically,  for  an  ADI  direction  of  length  N  we  need  only  12  parameters 
as  opposed  to  6N  parameters  for  nonuniform  meshes.  The  extra  6  parameters  are  for  t  e 
boundary  molecules.  A  further  advantage  is  that  for  a  given  problem  we  can  select 
coordinates  from  a  large  class  of  coordinate  systems.  In  the  process  of  sorting  throug 
the  various  possible  coordinate  systems  we  are  guided  by  two  criteria.  First,  the  new 
coordinates  must  lead  to  a  real  simplification;  and  secondly,  the  coordinates  must  be 
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easily  generated.  Since  bounding  surfaces  usually  become  coordinate  surfaces  the 
first  criterion  is  directly  measured  by  consideration  of  the  metric  tensor  (g^j)  which 
is  obtained  from  the  expression  for  the  fundamental  element  of  arc  length 


(ds)2  =  g  dy1  dyJ 


(1) 


Specifically,  an  increase  in  the  number  of  nontrivial  elements  in  the  expression  of  the 
metric  tensor  is  accompanied  by  a  corresponding  increase  in  the  number  of  terms  in  the 
equations  of  motion.  The  result  is  an  increase  in  the  computational  work  that  is 
needed  after  the  coordinates  have  been  generated  along  with  the  necessary  metric  data. 
The  second  criterion,  unlike  the  first,  is  most  often  neglected.  The  unfortunate 
result  is  that  there  is  often  more  work  involved  in  making  the  coordinates  than  in 
solving  the  original  problem  with  a  less  efficient  satisfaction  of  the  first  criterion. 
In  fact,  both  of  the  criteria  above  usually  are  at  opposite  polarities  in  complexity. 
The  prudent  selection  of  coordinates  is  then  a  balance  between  these  criteria. 


Our  criteria  for  selecting  a  suitable  system  of  coordinates  can  be  used  to 
compare  the  various  classes  of  coordinate  systems  and  to  evaluate  the  relative 
utility  of  each.  We  will  start  with  conformal  transformations  and  continually  enlarge 
the  class  until  we  obtain  general  none rthogonal  coordinates. 


For  conformal  transformations  the  metric  tensor  is  simply  given  by  a  scalar 
multiple  of  the  identity.  That  is,  g.,  =  My)  6^  where  the  kronecker  symbol  6^ 
vanishes  unless  i  =  j  in  which  case  itJis  unity.  From  this  expression  it  is  easy  to 

show  that  h  =  where  J  is  the  Jacob?  in  of  the  n-dimensional  conformal  trans¬ 

formation.  The  simplicity  of  the  metric  leads  to  very  simple  equations  of  motion  at 
the  expense  of  greatly  restricting  the  class  of  easily  obtained  transformations. 

These  transformations  are  generally  obtained  by  the  solution  of  partial  differential 
equations  which  may  in  itself  be  costly.  In  addition,  the  control  over  the  mesh 
distributions  is  indirect  at  best.  In  two  dimensions,  however,  conformal  trans¬ 
formations  have  been  successfully  used  on  many  occasions.  Here  the  metric  is  given 
by  g  =  J  j]  and  the  theory  of  functions  of  one  complex  variable  is  a  powerful 

tool  that  is  at  °  our  disposal.  When  the  boundaries  of  the  flow  region  can  be  matched 
with  well-known  conformal  transformations  there  is  nothing  that  can  effectively  compete 
with  this  way  of  generating  coordinates.  We  have  simply  optimized  on  the  generation 
of  the  equations  of  motion  and  on  their  solution  process  for  any  given  method  of 
solution.  In  a  number  of  cases  boundaries  can  be  matched  through  a  sequence  of  well- 
known  transformations.  However,  in  most  cases  of  practical  importance  the  boundaries 
are  too  complicated;  and  consequently,  cannot  be  simply  defined  as  desired. 


When  conformal  mappings  become  overly  difficult  to  construct,  it  is  best  to 
consider  the  slightly  larger  class  of  orthogonal  transformations.  For  orthogonal 
transformations  the  metric  tensor  is  given  by  the  diagonal  form  -  Pw)J  6^. 
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Note  that,  unlike  the  conformal  transformations,  the  diagonal  entries  of  the  metric 
can  be  different.  The  deviation  from  conformality  can  now  easily  be  measured  by  an 
examination  of  the  ratios  of  the  functions  To  see  why  this  is  so  we  need  an 
explicit  geometric  interpretation  of  the  metric.  For  a  position  vector  field  x,  the 
vector  field  e^  =  dx 
generated  by  holding 

It  is  often  common  practice  to  use  the  op^ator  notation  where  the  position  vector 
field  is  omitted.  By  an  application  of  the  chain  rule,  the  fundamental  element  of 
arc  length  can  be  expanded  as 


« 


iterpretation  of  the  metric.  For  a  position  vector  field  x,  the 
:/dp-  is  the  natural  tangent  vector  field  along  coordinate  curves 
the  remaining  coordinates  y  ,...,  y^""^">  y n  constant. 


(ds)2  =  dx.di?  =  (|4  =  ffl  *  = 


(2) 


and  hence,  by  linear  independence  g^j  =  e^.e. .  Now  note  that  we  have  an  orthogona- 
metric  if  and  only  if  the  e^  and  ej  are  perpendicular  when  i  f  J.  But  perpendicularity 
of  e*  and  ej  at  a  point  is  equivalent  to  the  perpendicular  crossing  of  the  associated 
coordinate  curves  at  the  point  in  question.  Consequently  our  intuitive  notion  of 
orthogonallity  in  terms  of  coordinate  curves  is  equivalent  to  the  metric  expression 
above.  In  addition  the  functions  are  easily  seen  to  be  equal  to  the  lengths  of 
the  corresponding  natural  tangent  vectors  % .  In  a  small  neighborhood  of  a  point 
the  functions ^  are  nearly  equal  to  their  values  at  the  point  and  thus,  the  measurement 
of  distance  along  coordinate  curves  is  very  nearly  given  by  distance  measurements  along 
the  respective  vectors  ei  in  the  tangent  plane  at  the  point  in  question.  When  the 
functions  tu  are  all  equal,  the  distance  measurement  in  the  tangent  plane  is  merely  a 
uniform  dilation  or  contraction  of  the  original  cartesian  system.  Thus,  length  ratios 
and,  hence,  angles  are  preserved  between  the  cartesian  system  and  the  tangent  plane. 

But ’the  projection  of  tangent  vectors  onto  the  curvilinear  system  preserves  angles. 
Hence,  with  equal  diagonal  entries  the  transformation  preserves  angles  and  is,  therefore 
called  conformal.  Consequently,  as  the  ratios  of  the  hj_  deviate  from  unity,  the 
transformation  smoothly  deviates  from  conformallity.  With  fewer  constraints  on  the 
metric  the  selection  of  coordinates  from  the  class  of  orthogonal  transformations  is 
slightly  less  restrictive  than  a  selection  from  the  class  of  conformal  transformations. 
The  process  of  coordinate  generation  is  usually  accomplished  by  geometric  methods  which 
result  in  a  system  of  differential  equations.  The  solution  of  these  equations  is  often 
too  costly  to  reasonably  perform  just  to  obtain  coordinates.  In  addition,  it  may  be 
difficult  or  even  impossible  to  properly  distribute  mesh  points  and  still  per  serve 

orthogonallity. 


General  nonorthogonal  coordinates  are  often  to  be  preferred  since  the  mesh  dis¬ 
tributions  can  be  controlled  and  since  the  coordinates  are  considerably  easier  to 
generate.  The  construction  process  is  entirely  geometric  and  generally  does  not  rely 
on  the  solution  of  differential  equations.  Furthermore,  points  can  be  essentially 
distributed  at  will.  Mesh  distribution  functions  can  often  be  directly  inserted  into^ 


the  transformation  in  a  manner  which  directly  distributes  the  points.  The  considerab 
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improvement  in  flexibility  associated  with  the  class  of  general  spacial  coordinates 
does  come  with  a  small  price.  Specifically,  the  metric  tensor  has  generally  non¬ 
trivial  off  diagonal  elements.  As  with  the  difference  between  orthogonal  and 
conformal  coordinates,  the  deviation  of  the  general  nonorthogonal  coordinates  from 
orthogonallity  can  be  measured  directly  from  the  metric.  That  is,  the  cosine  of  the 
angle  between  distinct  coordinate  curves  is  given  by  the  dct  product  of  the  associated 
unit  tangent  vectors.  The  cosine  of  the  angle  between  curves  i  and  j  can  be  written 
as: 


Thus  when  vanishes  for  distinct  i  and  J  we  have  orthogonallity,  and  when 
increases  from  o  the  coordinates  smoothly  deviate  from  orthogonallity  with  deviation 
given  by  the  arc  cosine  of  the  above.  This  deviation  can  be  used  to  advantage  by 
creating  almost  orthogonal  coordinates  in  certain  regions  of  importance.  For  example, 
one  may  wish  to  treat  boundary  layers  with  nearly  orthogonal  coordinates  and  let 
regions  of  greater  nonorthogonallity  fall  into  largely  inviscid  regions. 

I 

1 

With  all  of  the  above  considerations  born  in  mind  it  is  clear  that  the  general 
nonorthogonal  coordinates  are  the  most  suitable  choice  for  the  numerical  calculation 
of  the  fully  viscous  flow  field  over  an  ogival  body  at  angle  of  attack.  In  particular 
the  tube-like  coordinate  systems  to  be  discussed  in  a  later  section  are  general  non¬ 
orthogonal  coordinates  which  are  ideally  suited  to  the  ogival  body  problem.  Coordinate 
generation  with  tube-like  coordinates  is  computationally  efficient  and  at  the  same 
time  is  flexible  enough  to  allow  for  a  great  degree  of  control  over  the  mesh  distribution. | 
This  control  is  needed  for  the  resolution  of  boundary  layers  and  other  regions  where 
large  gradients  occur.  The  basic  choice  of  these  nonorthogonal  coordinates  essentially 
places  the  emphasis  on  the  fluid  dynamic  problem  rather  than  on  the  generation  of 
coordinate  systems  per  se. 


*  ■  i-i  -i  „  -  ~v.  ...  W..,  -- 
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THE  GENERAL  FORMATION  OF  AN  INITIAL  VALUE  PROBLEM  FOR  STRONGLY  CONVECTIVE  FLCWS 


Central  to  the  present  analysis  is  the  formulation  of  approximate  governing 
equations  which  can  be  solved  by  forward  marching  integration  in  the  direction  of  a 
"primary  flow".  The  entire  flow  field  can  thereby  be  obtained  by  a  sequence  of 
essentially  two-dimensional  calculations.  This  feature  of  the  method  results  in 
a  substantial  saving  of  computer  time  and  storage  compared  to  that  which  would  be 
required  for  solution  of  the  full  Navier-Stokes  equations.  The  equations  are  derived 
in  a  coordinate  independent  manner.  A  vector  field  that  reasonably  approximates  the 
primary  flow  direction  is  chosen  and  then  used  as  the  basis  for  an  approximation  of 
the  stress  tensor.  The  time-averaged  equations  are  written  in  general  conservation 
law  form,  and  then  the  approximate  stress  tensor  is  inserted  to  obtain  the  approximate 
equations .  Note  that  this  process  depends  only  on  the  choice  of  a  primary  vector 
field,  and  not  on  the  particular  coordinate  system  used  for  the  numerical  solution. 

The  primary  vector  field  used  here  consists  of  the  tangent  vectors  to  a  certain  family 
of  coordinate  curves  that  are  roughly  aligned  with  the  flow  geometry. 


The  governing  equations  are  derived  from  the  Navier-Stokes  equations  for  com¬ 
pressible  flow  of  a  viscous,  perfect  gas.  In  conservation  law  form  (Ref.  23)  and,  in 
general  curvilinear  coordinates  (y1,  y2,  y3),  these  equations  are  given  by 


H  t  A  (pvVe)  -  0 

a  t  a  y1 


for  continuity  and 


[  ovi  *4/7  1  +  — — »•  [(PvV+t1;J)  SxJ/g 

a  t  I  ay1  9yJ  By1 


for  momentum.  Constant  total  temperature  is  assumed,  and  thus  an  energy  equation  is 
not  required.  We  have  used  (x1,  x2,  x3)  for  fixed  cartesian  coordinates,  p,  for  density, 


v  =  v^  for  velocity,  g  =  det(gi^)  =  j  det  (—j)  j  for  the  metrical  determinant,  and 

for  the  components  of  the  stress  tensor  in  the  basis  e^  ®  e^ .  In  terms  of  the  metric, 
the  components  of  the  stress  tensor  are  given  by 


ij  ij  ,  ij  k  ,  o  ^ 

T  =  S  P  +  Q'k  v+ek  3e 


(5a) 
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where 


or 

k 

-•*<1 

(5b) 

and 

ild 

Pk 

-»<§ 

S1J  a* 

■  “  i  ■  ^  i  > 

(5c) 

for  viscosity 
symbols 

A  A 

|x,  inverse  metric  g  , 

i  a  ■* 

Kronecker  deltas  6j  =  6  d 

=  6^,  and  Christoff  el 

(5d) 


From  the  ideal  gas  law  and  the  constant  total,  temperature  assumption,  the  perfect  gas 
relation  has  the  form 


P  =  Ap  +  Bpg^vM 

where  A  and  B  are  constants.  In  all  of  the  above,  the  Einstein  summation  convention 
is  assumed.  That  is,  matching  upper  and  lower  indices  are  to  be  summed  from  1  to  3 
unless  otherwise  stated. 

It  is  assumed  that  for  high  Reynolds  number,  viscous  effects  are  negligible  except 
in  thin  layers  near  the  walls,  and  thus  boundary  layer  concepts  can  be  employed  to 
examine  the  relative  importance  of  viscous  terms  in  the  governing  equations.  Consequently, 
viscous  terms  which  are  considered  important  for  boundary  layer  flow  on  walls  are 
retained;  other  viscous  terms  are  neglected.  In  this  sense,  the  present  approach  can 
be  regarded  as  a  natural  extension  of  three-dimensional  boundary  layer  theory.  Unlike 
conventional  boundary  layer  theory,  however,  the  approximate  equations  are  to  be 
applicable  in  the  inviscid  flow  region  as  well  as  the  viscous  region  and,  thus,  no 
approximations  are  made  for  inviscid  terms  other  than  those  to  be  used  for  the 
pressure  field  in  subsonic  flow. 

To  account  for  turbulent  transport  processes,  the  governing  equations  are  time- 
averaged  in  the  usual  manner  for  turbulent  flows  (e.g.,  Hinze,  Ref.  24).  This  process 
of  averaging  produces  turbulent  correlations  which  are  conventionally  termed  Reynolds 
stresses.  Certain  components  of  viscous  stress  are  removed  from  the  time-averaged 
equations.  The  process  of  viscous  approximation  is  based  upon  the  assumption  that  a 
primary  flow  direction  exists.  This  direction  is  assumed  to  be  given  in  the  form  of 
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a  vector  field  ^  which  identifies  a  direction  at  each  spatial  point.  Then  can  be 
extended  to  an  orthogonal  triple  of  vectors  ?2,  ?3  a-fc  each  spatial  point. 

This  extension  must  be  accomplished  in  a  smooth  enough  fashion  over  the  whole  flow 
field  so  that  at  least  one  derivative  can  be  taken.  The  required  differentiability 
occurs  because  of  the  requirement  to  differentiate  the  components  of  the  stress  tensor 
as  they  occur  in  the  Navier-Stokes  equations.  From  this  construction  one  obtains  a 
specification  of  a  field  of  orthogonal  frames  where  ^  is  aligned  with  the  assumed 
primary  flow  direction,  and  for  each  point  and  ?2  span  an  orthogonal  transverse 
plane.  Within  such  a  frame,  we  form  a  differential  viscous  stress  cube  (Fig.  2). 

The  resulting  components  of  viscous  stress  on  the  cube  surface  are  either  aligned 
with  or  are  orthogonal  to  the  primary  direction.  In  this  way  the  force  balances 
represented  by  the  Navier-Stokes  equations  are  effectively  separated  into  three 
mutually  exclusive  directions  so  that  approximations  in  any  given  direction  do  not 
directly  affect  other  directions.  That  is,  forces  in  any  one  of  the  directions  do 
not  have  nontrivial  projections  on  the  other  two  remaining  directions.  If  the 
equations  of  motion  were  written  for  the  isolated  cube,  then  the  stress  components 
would  contribute  to  the  force  balance  in  their  respective  directions.  However, 
in  the  primary  direction  the  viscous  contribution  o*33  is  expected  to  add  little 
to  the  strong  convective  forces  and,  hence,  this  contribution  is  ignored.  In  addition, 
the  contribution  of  the  viscous  shearing  stresses  a31  and  o3  to  the  force  balance 
in  the  transverse  equations  for  the  5^  and  ?2  directions  is  also  small  relative  to 
convective  forces.  These  force  balances  are  mutually  exclusive  due  to  the  orthogonality 
of  the  frame.  Effectively,  we  can  generate  longer  and  longer  viscous  stress  cubes  by 
joining  existing  cubes  along  transverse  faces.  The  total  assumption  is  that  viscous 
forces  on  transverse  faces  are  negligible.  Eventually,  we  are  considering  forces  on 
a  fiber-like  object  aligned  with  the  primary  direction  (Fig.  2).  From  this  viewpoint, 
we  are  just  neglecting  internal  viscous  forces  within  the  fiber.  That  is,  the  fiber 
has  no  stiffness  and,  therefore,  the  only  balance  against  the  convective  forces  is  due 
to  the  shearing  stress  along  its  boundary.  This  is  particularly  appropriate  when  the 
fiber  is  in  the  boundary  layer,  where  the  no-slip  condition  causes  the  fluid  to 
decelerate  and  come  to  rest  at  the  walls. 

For  a  viscous  stress  cube  in  the  frame  £2,  §3,  the  viscous  stresses  have  the 

form 


(6) 


and  it  is  postulated  that  the  components  o  for  i  =  1,2,3  are  negligible  relative  to 
convective  forces.  Thus  we  shall  replace  the  tensor  a  with  the  tensor 

g)  Fj  where  a  =  (1-63)  o^.  Unlike  the  original  tensor,  this  new  tensor  is 
not  symmetric.  When  this  is  inserted  into  the  time-averaged  Navier-Stokes  equations, 
we  obtain  an  approximate  set  of  governing  equations  which  are  not  elliptic,  and  which 
can  be  solved  by  a  forward  marching  procedure. 
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The  approximate  governing  equations  are  obtained  in  terms  of  the  metric  data 
associated  with  a  curvilinear  coordinate  system.  In  the  ogival  body  problem,  tube¬ 
like  coordinates  are  used  to  match  coordiate  surfaces  with  the  ogival  body  and  the 
bow  shock.  The  solution  is  obtained  by  a  forward  marching  procedure  down  the  axis 
of  the  ogival  body.  The  development  and  analysis  of  tube-like  coordinates  is  presented 
in  a  later  section,  together  with  the  metric  data.  For  the  present,  we  will  assume 
that  we  are  given  an  arbitrary  coordinate  direction  taken  as  a  marching  direction  for 
the  solution.  This  direction  need  not  be  aligned  with  the  primary  flow  direction  5^. 
The  natural  basis  of  vectors  tangent  to  coordinate  curves  is  given  by  e±  =  d/d y  for 
i  =  1,2,3.  But  the  basis  |lf  ?2,  for  our  approximations  is  independent  of  this 
coordinate  basis.  Hence,  our  approximation  scheme  is  executed  by  transformation  of 
the  viscous  stress  from  the  e^  basis  into  followed  by_the  stress  approximation  as 
described  above  and,  finally ,_a  transformation  from  the  Ij.  basis  back  into  e^  This 
requires  the  transformation  £.  =  and  its  inverse  e^=  Navier-Stokes 

equatios,  the  components  of  stress'5  are  represented  by  a  sum 


Tij  =  pg1^  + 


(7) 


where  is  the  viscous  part  of  the  stress  tensor.  But  the  r5-*5  are  just  the 
coefficients  of  the  stress  tensor  in  the  natural  basis  of  tangents  to 

Thus  we  must  express  ius  viscous  part  in  terms  of  the  frame 


coordinate  curves. 
-* 

?1  >  ^>2*  *3* 


S  ,  So*  By  a  change  of  tensor  product  basis  we  obtain 


viJ 


ej  = 


ij 


x  a 
\ 


L©  5 


(8) 


Now  the  coefficients  in  the  tensor  product  basis  of  frame  vectors  are  replaced  by 
their  approximates.  This  yields 


-  0-»?)  Til  1?  ?r®rs 

-  A $  15  r®?s 


(v«  Ti\  Hj)  ?r®?s 


(9) 
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A  transformation  back  into  the  tensor  product  basis  of  the  coordinate  tangents  is 
needed  so  that  the  approximation  can  be  properly  inserted  into  Navier-Stokes  equations 
which  were  expressed  in  terms  of  the  coordinate  system.  We  obtain 


(»1J  <lfi)  5© 5 

-L  J 


=  £  T|r  T|*  5*  §s  ©  4 

r=l 


4 i  r  „m  „k 

*  rEi  « J  Hi  5r  ^®«k 


2 

£ 

r=l 


V  "Lie  -*  -♦ 

-  v  \  em®  ek 


(10) 


and,  therefore,  new  coefficients  of  viscous  stress  u>  defined  by 


mk  ~  mk 
0)  =  v 


ik 

v 


(11) 


-3  -m  >  <  ,r  J  .  n  ^  ^ 


<  m  _3  m  ,  <  ik  j 

=  W  -  h  1  I  “i  v  +  *< 


I 


•i  "1  ’3  >  '  j  -3 

The  approximation  is  now  complete  when  this  is  inserted  into  the  equations. 


As  a  special  case,  suppose  that  the  primary  direction  vector  field  is  given  by 
the  y3  coordinate  curves.  Then,  =  ®3*  This  is  eaFily  extended  to  an  orthogonal 
triad  by  setting  |g  =  e2  and  =  g  e  In  the  present  tube-like  geometry,  we 

■i  2  q  31  3 

would  take  y  ,  y  ,  and  y*5  for  the  angular,  radial,  and  axial  variables,  respectively, 
as  illustrated  in  Fig.  3. 


Our  approximation  to  the  Navier-Stokes  equations  can  be  written  explicitly  in 
the  conservation  law  form  given  previously.  The  time-averaging  for  turbulent 
fluctuations,  however,  requires  some  additional  notation.  Specifically,  the  dependent 
variables  are  represented  as  the  sum  of  a  time-averaged  quantity  denoted  by  a  straight 
overbar  (-)  and  an  instantaneous  fluctuating  quantity  denoted  by  a  curved  overbar  (~) . 
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After  time-averaging,  the  equations  become 


Cp  v*  +  T  v1  )  /g 


=  0 


(12a) 


for  continuity,  and 


+  p  yi  VJ  +  p  +  -p 

+  gi^  p  +  »«)  ~^[  /g|  = 


v1 


0 


(12b) 


for  momentum.  The  triple  cross  correlation  terms  have  been  neglected.  Since  the 
remaining  components  of  Reynolds  stress  are  coefficients  in  the  tensor  product  basis 
e.  ®  e  ,  they  can  be  expressed  in  any  needed  form  via  a  change  of  tensor  product  basis. 

Note  that  since  the  new  stress  tensor  is  not  symmetric  the  index  ordering  in  the 
viscous  stress  is  important.  In  fact,  if  by  mistake  it  should  become  inverted,  then 
the  primary  direction  momentum  balance  would  have  no  viscous  contributions  at  all. 
This,  however,  would  be  a  false  result  since  significant  shearing  stresses  would  be 
neglected.  One  can  most  easily  examine  this  situation  via  the  stress  cube  (Fig.  2) 
by  considering  any  cartesian  system  with  one  coordinate  aligned  with  the  primary 
direction  With  proper  index  alignment  one  gets  precisely  the  derired  equations. 

However,  with  the  incorrect  alignment  the  axial  momentum  equation  is  inviscid.  Thus 
considerable  care  must  be  taken  so  that  one  does  not  inadvertantly  apply  symmetry  to 
the  nonsymmetric  stress  tensor  derived  here. 

For  entirely  supersonic  flows,  the  approximate  equations  (Eqs.  (12a)  and  (12b)), 
together  with  boundary  and  initial  conditions  can  be  solved  by  forward  marching 
integration  in  the  x  direction  without  any  assumptions  about  the  pressure  field,  as 
was  demonstrated  by  McDonald  &  Briley  (Ref.  18)  for  laminar  flow  in  rectangular  jets. 
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THE  NUMERICAL  METHOD 


Overview 

In  the  previous  section  we  derived  a  system  of  equations  in  which  the 
time-like  derivative  of  the  solution  vector  was  expressed  implictly.  This  occurred 
because  time-like  derivative  was  applied  to  a  flux  vector  instead  of  the  solution 
vector  directly.  The  association  with  the  word  flux  is  a  direct  result  of  the 
general  Stokes  theorem  and  the  conservation  law  form  of  Eqs.  (12a)  and  (12b). 

Since  the  time-like  direction  is  the  spatial  marching  direction  and  since  the  equa¬ 
tions  were  derived  from  the  steady-state  Navier-Stokes  equations,  (4a)  and  (4b), 
the  flux  vector  is  a  nontrivial  function  of  the  solution.  Consequently  one  is 
faced  with  the  problem  of  constructing  a  numerical  method  that  is  sufficiently 
general  to  solve  the  general  system  of  equations  of  the  previous  section.  Since 
time-like  variations  are  nontrivial,  a  second  order  solution  is  preferred. 

The  method  used  is  based  on  an  implicit  scheme  which  is  potentially  stable  for 
large  step  sizes.  Thus  as  a  practical  matter,  stability  restrictions  which  limit 
the  axial  step  size  relative  to  the  transverse  mesh  spacing  and  which  become  pro¬ 
hibitive  for  even  locally  refined  meshes  (e.g.,  in  laminar  sublayers)  are  not  a 
factor  in  making  the  calculations.  The  general  approach  is  to  employ  an  implicit 
difference  formulation  and  to  linearize  the  implicit  equations  by  expansion  about 
the  solution  at  the  most  recent  axial  location.  Terms  in  the  difference  equations 
are  then  grouped  by  coordinate  direction  and  one  of  the  available  alternating- 
direction  implicit  (ADI)  or  splitting  techniques  is  used  to  reduce  the  multidimen¬ 
sional  difference  equations  to  a  sequence  of  one -dimensional  equations.  These 
linear  one-dimensional  difference  equations  can  be  written  in  block- tridiagonal  or 
a  closely  related  matrix  form  and  solved  efficiently  without  iteration  by  stan¬ 
dard  block  elimination  techniques.  The  general  solution  procedure  is  quite  flexible 
in  matters  of  detail  such  as  the  type  and  order  accuracy  of  the  difference  approxi¬ 
mations  and  the  particular  scheme  for  splitting  multidimensional  difference  approxi¬ 
mations.  Based  on  previous  experience  of  the  authors,  however,  it  is  believed  that 
the  consistent  use  of  a  formal  linearization  procedure,  which  incidentally  requires 
the  solution  of  coupled  difference  equations  in  most  instances,  is  a  major  factor 
in  realizing  the  potential  favorable  stability  properties  generally  attributed  to 
implicit  difference  schemes. 

After  modeling  Reynolds  stresses  the  governing  equations  of  the  previous 
section  can  be  solved  by  the  general  approach  to  be  presented  herein  which  is 
motivated  by  the  previous  work  of  McDonald  and  Briley  (Ref.  25).  Unlike  the  earlier 
work,  the  present  algorithm  achieves  second  order  accuracy  for  the  forward  marching 
direction  using  only  two  levels  of  storage  as  opposed  to  three.  In  addition, 
the  present  analysis  is  developed  directly  on  tensor-like  objects  which  contain 


17 


R76-912024-8 


a  considerable  amount  of  generality  and  hence  applicability.  Further  developments 
come  with  the  construction  of  an  implicit- explicit  shifting  function  which  can  be 
used  to  well-condition  an  ill-conditioned  transformation  by  creating  a  diagonally 
dominant  (Varga,  Ref.  26)  and  hence  numerically  solvable  system.  Other  applications 
of  implicit-explicit  shifting  occur  when  it  is  desired  to  cast  the  implict  differ¬ 
ence  equations  in  a  form  suitable  for  ADI  splitting  along  coordinate  directions. 
Here  one  shifts  all  mixed  implict  transverse  derivatives  to  the  explicit  side. 

This  avoids  extra  ADI  sweeps  or  an  increase  in  the  matrix  banded  structures  of  the 
individual  sweeps. 


The  System  of  Equations 


In  keeping  with  the  viewpoint  of  maintaining  generality,  the  numerical  method 
shall  be  developed  in  a  slightly  more  general  context  than  is  actually  needed. 
Specifically,  M-th  order  N-dimensional  systems  of  the  form 


6% 

6t 


F< 


(13a) 


are  to  be  solved  for  a  solution  (cp^)  =  (<p]_,  .  .  .  cpjj)  where 


Hi  =  H^t,  Xj,  Fi0f)  =  Hi(T,  Xj,  Picr(X,  xk,  tPjg)) 
Fi  =  F i(r,  Xy  qia)  =  Fi(T,  Xy  qiaU,  xk,  cp^) ) 
T(t)  =  X ( t)  =  t  =  time -like  variable 
(Xj)  =  (XL,  .  .  .,  Xjj)  =  spatial  variables 


(13b) 

(13c) 

(13d) 

(I3e) 


for  l£i,  j,  k,  and  a  multi-index  a  =  (<*]_,  .  •  •>  «2m)  ln,fceSers  “j  • 

The  reader  should  note  that  there  has  been  a  use  of  a  compressed  notation  to 
indicate  functional  dependence.  The  power  of  this  notation  is  easily  observed  from 
the  expansion. 


piQf(X,  xk,  cp*)  =  PiQ?(X,  x1#  xg,  .  .  .,  xN,  <pL,  <p2,  •  •  •>  'Pn) 


(14) 


In  addition,  this  notation  avoids  the  confusion  that  can  result  when  chain  rule 
expansions  produce  Jacobian-like  objects  of  varying  dimensions.  This  should  be 
obvious  from  the  form  H.(r,  Xy  Pia)  since  (xj)  and  (picr)  are  vectors  of  generally 
unequal  lengths.  This  notation  should,  of  course,  be  distinguish^  from  the 
Einstein  summation  convention  of  summing  like  indices  which  is  iu  the  context  of 
sums  of  products  as  opposed  to  the  argument  listings  here.  The  spatial 
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derivatives  and  the  dependence  upon  the  solution  vector  (cp^)  are  all  contained 
in  the  specification  of  the  functions  pia  and  qix  in  (I3h)  and  (13c) .  Specifically, 
one  has 


where 


and 


p,  =  D  G  .  •  •  • 

xa 


G. 


xa, 


2M 


(15a) 


Identity 


for  a.  *  k,  2, 

for  a.  =  0 
J 


(15b) 


G, 


lQto+j 


=  G.  (X, 
laM+j 


ck’  fy) 


(15c) 


The  form  of  qiff  is  similar  but  with  different  functions  Gj^  .  It  should  be 
noticed  that  all  of  these  specific  forms  can  be  expanded  by  successive  applications 
of  the  chain  and  Leibnitz  rules.  If  these  expansions  were  taken  to  the  point  where 
only  derivatives  of  the  form  wou^<i  aPPear  then  both  and  qiQ,  could 

be  replaced  by  expressions  of  the  form  D^cp^  where  1^  =  1^  . .  .I^M  and  a  =*  (o^, 

Qf„,  a  , ,  0,  0)  with  0  so.  sn  and  ,  f  0.  In  particular,  one  would  have 

the  solution  vector  3  Pia  3  for  »  3  (0, . . .  ,0,a  0,...,  0)  and  deriva¬ 
tives  of  uP  bo  order  M  as  a  varied  otherwise.  This  simplified  functional  fora 

for  pia  and  q^  would  certainly  make  the  theoretical  discussion  easier;  however, 

the  computational  complexity  of  the  problem  would  generally  be  increased.  In 
many  cases  a  natural  grouping  of  functions  p^a  and  q^a  occurs.  Such  groupings  are 
easy  to  spot  in  the  general  fluid  equations  (12a)  to  (12b).  The  above  expansion 
and  redefinition  of  p^a  and  q^a  would  generally  increase  the  number  of  terms  in  the 
equations .  Each  of  these  terms  would  require  roughly  the  same  number  of  operations 
involved  in  linearizations  and  differencing  as  was  required  in  the  original.  Thus, 
the  operation  count  would  generally  go  up  in  a  direct  proportion  to  the  increase 
in  number  of  terms .  Consequently,  it  is  preferrable  to  stay  with  the  theoretically 
more  cumbersome  form  of  equations  (15)  which  produce  a  computationally  more  efficient 
scheme.  Since  the  derivative  operators  with  respect  to  X,  <Pg,  and  t  all  commute 
with  the  operator  .  of  (15b);  the  differentiation  of  p^  and  q^,  with  respect  to 

X,  9^,  or  t  is  easily  seen  to  obey  the  chain  rule.  Thus,  one  can  treat  the  operator 

functions  p^a  and  q^  as  ordinary  functions  in  the  variables  X  and  cp^. 
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From  the  chain  rule,  the  system  of  partial  differential  equations  (13)  can  be 
rewritten  in  the  form 


A  =  F  _  ^  — -  (l6a) 

i£  ftt  1  &t  bPia  bX 


where 

.  B  =  SL  (l6b) 

bPia  b'Pjj 


is  just  the  Jacobian  of  the  flux  (Hi)  in  the  solution  vector  (<PA).  If  the  matrix 
(A  )  represents  a  nonsingular  linear  transformation  A,  then  one  can  directly  solve 
for  the  time-like  derivatives  of  the  entire  solution  vector.  Otherwise  one  must 
consider  a  system  of  lower  rank  and  with  constraining  relations.  Under  a  change  of 
basis  equations  (l6)  can  be  rewritten  in  an  equivalent  form  where  the  linear  trans¬ 
formation  A  is  represented  by  a  matrix  with  r  linearly  independent  rows  and  N-r 
rows  of  zeros  where  r  is  the  rank  of  A.  The  last  N-r  rows  of  the  transformed  equa¬ 
tions  correspond  to  the  nullity  of  A  and  hence  involve  no  time-like  derivatives. 
Consequently  one  can  consider  these  to  be  constraining  relationships  which  can  be 
used  to  eliminate  the  last  N-r  rows  of  the  solution  vector  in  the  first  r  equations. 
If  this  is  possible,  then  one  has  reduced  the  original  system  of  N  equations  in 
N  unknowns  to  an  equivalent  system  of  r  equations  in  r  unknowns.  These  systems  will 
be  called  reducible  systems.  This  terminology  corresponds  with  the  matrix  terminol¬ 
ogy  since  a  discretization  of  the  system  would  lead  to  reducible  matrices  (Varga, 
Ref.  26).  If  the  system  is  reducible  at  each  point  (t,  x. )  with  a  continuous  trans¬ 
formation  of  basis,  then  the  system  will  be  called  a  solvable  system.  Only  solvable 
systems  shall  be  considered  here;  and,  in  fact,  without  loss  of  generality  one  can 
assume  that  A  is  nonsingular  for  otherwise  the  reduced  system  would  just  produce  a 
smaller  nonsingular  version  of  A  which  would  be  solved  in  the  same  manner.  With 
this  consideration  one  obtains 


bHj  _  bHi__  bPio,  | 

?>T  bPia  1 


(17) 


where  the  matrix  (B-h)  is  the  inverse  of  (Aik)  and  the  notation  Sj  is  introduce 
to  denote  the  value  of  the  time-like  derivative  of  the  jtn  component  of  the  solu¬ 
tion  vector  that  is  determined  directly  by  the  system  of  partial  differential  equa¬ 
tions  . 
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The  Numerical  Scheme 

Time-like  variations  of  the  system  of  equations  (13)  are  generally  nontrivial. 
This  is  evident  on  examination  of  the  system  of  approximate  Navier-Stokes  equations 
developed  in  the  previous  chapter.  There  one  has  a  time-like  direction 
associated  with  the  spatial  direction  which  is  then  used  to  march  the  solution. 

The  remaining  directions  form  what  is  generally  referred  to  as  a  transverse  surface 
which  is  often  taken  as  a  transverse  plane.  As  the  solution  is  marched  the  geo¬ 
metry  of  the  problem  changes  as  the  bounding  contours  in  transverse  surfaces 
generate  the  flow  region  in  a  generally  nontrivial  fashion.  Thus  one  has  con¬ 
siderably  more  opportunities  for  error  growth  than  if  the  flow  region  had  no  varia¬ 
tions  in  the  time-like  direction  at  all.  With  this  motivation  it  is  preferrable 
to  develop  a  second  order  scheme  as  opposed  to  a  first  order  one. 

The  numerical  scheme  that  is  developed  here  is  a  second  order  generalization 
of  the  classical  Crank-Nicolson  scheme.  In  the  well-centered  framework  of 
Crank-Nicolson  one  has 


n+2  / 

=  p  +  O(h^) 

i 


where  h  =  tn+1  -  tn  and  all  superscripts  are  used  to  denote  time-like  evaluations. 
The  reader  may  also  notice  that  all  other  indices  in  this  section  were  carefully 
taken  to  be  subscripts.  In  this  way  no  confusion  can  result.  One  now  must  evalu¬ 
ate  both  sides  of  the  Crank-Nicolson  equation  (l8)  and  preserve  second  order 
accuracy.  From  a  Taylor  expansion  of  the  right  hand  side  about  level  n  one  obtains 


n+^  -  F  n  +i££i  +  +  1  f  £ 

1  1  l  &T  bqia  1  bX  b<Pj>  bt  M  2 


+  0(h2) 


(19) 


where  the  chain  rule  has  been  used.  The  n-level  evaluations  in  the  first  order 
piece  are  straight  forward  with  the  exception  of  the  quantity  (bP^/bt)11.  This  can 
be  either  evaluated  by  a  finite  difference  or  directly  from  the  differential 
equations  with  Sjn  of  equation  (17)  •  IT  the  latter  approach  is  taken,  then  the 
implicit  character  of  the  basic  Crank-Nicolson  scheme  is  lost.  Thus  a  finite 
difference  shall  be  used.  Since  the  term  itself  is  first  order,  the  simple  first 
order  forward  difference  is  sufficient.  Thus  one  obtains 
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*Fi  J  [b^}  *  5 


'bT  '  'Mia  0K 

MT.  xn  ,xn.  \n  n+1 

wHi_.  £5ia  (vt 

'•V  'tffi  1 


nv  /  2. 

)  +  0(h  ) 


On  the  left  hand  side  of  the  Crank-Nicolson  equation  (l8)  one  must  evaluated  a 
time-like  difference  quotient  of  fluxes.  This  is  obtained  by  a  Taylor  Series 
expansion  about  the  n-level.  To  break  the  calculation  up  into  manageable  pieces 
it  is  best  to  first  expand  each  of  the  (n+Jr)- level  derivatives  about  level  n.  I 
so  doing  one  obtains  the  finite  difference  « 


(£) 

the  Taylor  Series  expansion 


iL.  +  0(h2), 


tv^F1 , 

(*k)\  ( 

'bT  ' 

bT  ' 

-(£>  *  S1* 


2  ' 

j  &^i  +  &2Hj 
I  6t2  bPi(yJ)i 


n 

5Pi2+  h 

b\  b<Pi  bt  2 


+  0(h2) 


and  the  similar  expansions 


/3l) 

Ld,  'hP.  ' 


bTbV.  ‘ 
ia 


i _  [i£i£  +  ^i|  [n  | 

.ftp.  bX  b<P,  bt  J  > 


bV.gbV. 
i  ia 


+  0(h2) 


L,  '  y  b\  ' 


ha  [  &  pjq 
-  b<PjbX  bt 


|  +  0(h2) 


and 
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n+i  ^ 

(Sat)  -(2£te)  + 

>*  7  '*»* 


b2p 


2  n 
b  Pia 


bXb^  b'Pjb^ 


bt 


I  +  0(h2) 


(21e) 


The  time-like  difference  of  equatior  (18)  is  now  obtained  by  a  direct  substitution 
of  equations  (21)  followed  by  an  absorbtion  of  second  order  terms  into  the  collec¬ 
tion  0(h2).  The  result  is  given  by 


n+1  n 

H.  -H. 

n 

_(B>) 

+  S^Mi 

+  b2Hi 

h 

Ibr2 

bPiabT 

*(—) 

'  bP- 
utia 

l'f*> 

n  x2_ 

+  1  LEifl 

l  bX2 

.in* 

b'P.bX 

M 

+  ( +  1 
Vb^  1 

2 

i  b  Pia 

1  bX&fy 

2 

,  b  Pic* 

b'P.b'P. 

3  l 

SH*) 

+jb2%  + 

b2}h 

bPi/5  b«Pj  | 

(  bTbP. 

Iff 

bPiebPia  ‘  bx 

b^  bt 
«J 

n  n-t  j. 

n 

\ 

,  /  &pio 

\  t*i 

-  ^  \  + . 

)  \  h 

Is 

n 


.  \Eis  + 

Wia'  *>.1 

1  bX 

WPj  bt 

n 

t  2ii 

h 

.  bt  1 

2 

n  n+1 

n 

1  V 

4 ) 

)  V 

h  M 

n 

n 

i!  -  ( 

bPj»\ 

li  2  H 

bX  ' 

(22) 


Now  it  is  best  to  regroup  the  above  into  n-level  time-like  derivatives  of  9.  and 
time-like  differences  in  <P£  with  coefficients  ordered  by  powers  of  h.  This  regrouping 

yields  the  foma 


Hi 


n+1  n 
-  Hi 


n+1 


n 


h 


n  „/b'Mn  n,<Pt  -  9*1  »  (£El f(fi 

‘i  '+  »>1J  (-sf)  +  CU  1 - Z  >+  dW*  lb*  >  1 


n+1  n 
"  *±) 


h 


+  0(h2) 


(23) 


where 
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and 


_^i  +  Wi  ftpiff  „  h  |6Hi  +2 

1  bT  bp.^bX  2  16t2 


b2Hj 


&Pia6r  &x 


bp^  +  &Hj 
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b  P 


i2 


&Pia  &X" 


b2Hi 


6piB*i«  bX 


SPiB  dPicy 

ax  ■ 


ij 


[cK.  b2pi(y 

-  ( 

b2Ht 

2 

b  Hi 
■f 

*ial 

lsPi(y  SXb<Pj 

bTbPicy 

bP. 

iB 

bP. 

iff 

bX  / 

1 

C  = 

bHi 

+ 

b, „ 

ije 

bPiu 

ije 

fbHi  b2p 

i a  + 

2 

b  Hi 

_®i£ 

bPfff 

u 

l  bP^Q,  bcpj 

bcp^ 

bpi0bPia  bcp^ 

j 2 

(24a) 


(24b) 


(24c) 


(24d) 


From  successive  applications  of  the  chain  rule  in  conjunction  with  the  Leibnitz  rule, 
the  coefficients  can  be  condensed  into  the  simple  forms 


bHi  h  b2Hi 

(25t% 

1  bT  ^  bT2 

h  b%i 

(25b) 

bij  2  bTb<Pj 

*'bHi 

c j  4  =  - +  b.  . 

iJ  b^j  iO 

(25c) 

2 

h  b  Hf 

ijk  ^  &cPjbVk 

(25d) 

where 

b_  =  b_  +  b_  * 
bT  ~  bT  bX 


(25e) 


Time-like  derivatives  of  quantities  other  than  the  solution  vector  cp^  are  given 
by  derivatives  with  respect  to  T.  Such  quantities,  for  example,  are  often  items 
that  are  built  up  from  the  metric  information  when  the  equations  of  motion  are 
expressed  in  some  curvilinear  coordinate  system.  The  derivatives  with  respect  to 
T  occur  only  in  the  coefficients  ay  by,  and  Cy.  When  the  time-^ike  dependence 
is  on  the  solution  vector  cp.  alone,  the  coefficients  ai  and  by  both  obviously  vanish 
and  the  coefficient  Cy  is  just  the  Jacobian  of  Hy  The  expression  for  dyk  is 
independent  of  T -derivatives  and  is  easily  identified  as  the  product  of  h/2  and  the 
Hessian  of  Hy 


t? 


f 
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m 

It; 
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To  maintain  second  order  accuracy  the  derivative 


5Sl 

bt 


associated  with  the 


d.  .  -terms  in  equation  (23)  cannot  be  evaluated  with  a  forward  difference;  for 
then,  the  quantities  V^n+1  would  appear  in  products  of  the  form  cpjn+1cpjJn+1  and  thus 
linearity  would  be  lost.  Consequently,  one  is  left  with  two  options.  First,  one 
may  take  a  simple  backwards  difference  to  a  level  (n-l)  in  the  time-like  variable. 


This  difference  represents  the  derivative  to  first  order;  and  since  the  coefficient 
d— .  already  contains  a  factor  of  h,  one  has  obtained  second  order  accuracy.  This 
was  the  option  taken  by  McDonald  and  Briley  (Ref.  25).  rm~~  1  '"~1" 

however,  a  three  level  scheme. 


The  resulting  scheme  is, 

In  some  problems  storage  is  the  critical  factor; 
and  in  such  cases, a  two  level  scheme  is  much  more  desirable  than  a  three  level  one. 
This  motivates  the  second  option.  The  idea  is  simply  to  evaluate  the  derivative 
directly  from  the  orignal  system  of  differential  equations.  Thus  from  equation 
(17)  one  has 


bCH 


Sj  bt 


=  B 


ji 


Fi  ~  bT 


(26) 


in  terms  of  the  T-derivative.  From  a  direct  substitution  of  this  expression  into 
the  d^j^- terms  of  equation  (23)  one  obtains  the  representation 


n+1  n 
H,  -  H 


i  "  Hi  n  ,  n  /b<Pj\  n  n 

h  i  iJ  \ht  /  if 


n+1 


n 


bt 


+  d.  .  s 
ijf  j 


n  „  n)  -V  j 


(27) 


+  0(h2) 


Under  either  of  the  above 


for  the  time-like  flux  difference  on  the  left-1  \nd  side, 
options  one  has  not  yet  finished  the  representation  of  this  time-like  flux  difference 
since  the  n-level  derivatives  )n  associated  with  the  b-11 

bt  n 

to  be  considered.  Since  coefficients  b,.  *  contain  a  factor  of  h  one  can  represent 


-terms  still  remain 


Since  coefficients  b^ 
these  derivatives  by  either  a  first  order  difference  or  a  direct  substitution  of 
sr.  Unlike  the  d^  .-terms  there  are  no  (n+l) -level  complications  here  and  conse 


j  |  U11X4.4LV 

quently  only  the  simple  first  order  forward  difference  need  by  considered.  The 

For  any  n-level  coefficient 


representation  of  these  terms,  however,  is  not  unique. 


e^j  that  is  in  0(h)  one  has 


'13 


0  + 


o(h2) 


(28) 


which  can  be  added  to  equation  (23)  without  changing  the  order.  Now  if  one  uses 
the  forward  difference  for  the  b^j- terms  and  adds  equation  (28)  to  equation  (27), 
then  one  obtains  the  time-like  flux  difference  representation 


'niiiir  iif  i  ii  iiTMMiaitiiMiir  m  iriniii 
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n+1  n 

H.  -  H. 


=  a.  -  e. 
i  i  i 


n+l  n 

i  i  I  i  n  n  ,  i  °1  ft  -  t  \ 

s  +  b .  +e.  +  c  \  / 

\  5 1  l  IjJ  ijl  0  M  Yi  ' 


+  0(h2) 


With  e.  n  =  0  the  effect  is  to  represent  b^-terms  with  a  finite  difference;  and  with 

O'  Y\  .  ^  1J _  J_T _ r.  v.n  rj  V»  A  Y\  PA  Tl  A 


e  n  -  -b-  n  the  direct  substitution  of  s/1  is  used.  In  between  there  are  an  infinite 

C*  4  —  U  T  ft  "  .  .  .  «  II _ _ _ •  J  ^  T.fi  4-V»  Q 


0,  —  ;=  —Q  .  .  OllC  J.  ~ - -  y-j  .  ,  , 

number  o/possibilities.  In  effect  this  arbitrariness  of  eijt  provides  one  with  a 
degree  of  freedom  that  can  be  used  to  shift  terms  between  implicit  and  explicit  parts 
of  tie  representation.  The  advantage  here  is  that  terns  can  be  shifted  towards  the 
implicit  side  to  insure  diagonal  dominance  and  hence  invertibility  of  the  implicit 
scheme.  In  addition,  unwanted  implicit  terms  can  be  shifted  to  made  ADI-splitting 
easier!  For  these  rLsons,  the  function  ei£n  shall  be  called  an  implicit-explicit 

shifter . 


Both  sides  of  the  Crnnk-Wicolson  equations  (18)  hare  been  suitably 
to  second  order  with  right  and  left  hand  sides  being  given  by  equations  (SO)  and 
(29)  respectively.  By  direct  substitution  one  has  the  second^order  scheme 

n  n  n  r  n  n  n  ,  n  nl  _ 


n  n 

eU  s*  + 


+  e. 
i  SL 


where  *  n+1  =  cp  n+1-cp/\  From  a  multiplication  by  h,  the  notation  of  equation 
(25e) ,  and  the  chain  rule  one  obtains  the  form 


*ii  +  eU  +  cu 


h  BUI*1  n+1 

+  dij£  SJ  "  2  6cp  j  h 


f  h  6Fi  In 

=  hlFi+i^_ai"e^  S^l 


where  the  strength  of  the  implicit  coupling  is  controlled  by  the  first  order  f^nc^±on 
e  .  While  this  function  can  be  taker  to  extreme  settings,  it  is  usually  taken  as 
zero  unless  mixed  derivatives  occur  or  diagonal  dominance  is  lost.  The  spatial 
differences  are  easily  inserted  into  equation  (3D  by  the  mere  replacement  of  the 
spatial  differential  operators  Dk  of  equation  (15b)  with  the  difference  operators 
which  are  defined  by 


Ak  =  /  1="r 


E  Wk(i)  Ek(D  for  o  <  k  < 


Ek(o) 


for  k  =  o 
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where  Ek(i)  is  the  evaluation  map  on  the  position  designated  by  translation  of  the 
kth  grid  variable  by  i  grid  points  and  where  Wk(i)  is  the  associated  finite  difference 
weight.  One  easily  notes  that  Ek(o)  is  the  identity.  When  the  weights  Wk(i)  are 
taken  to  give  second  order  spatial  accuracy  in  the  interior  and  on  the  boundary  the 
overall  scheme  of  equation  (31)  is  second  order.  Second  order  accuracy  at  interior 
points  with  a  uniform  mesh  is  achieved  by  setting  r=s=l,  Wk(o)=o  and 


Wk(l)  =  -  Wk(-1)  = 


(Ax). 


(33) 


while  taking  Ek(+l)  to  be  half  point  evaluations  via  two  point  averages, 


The  scheme  of  equation  (31)  produces  a  linear  implicit  system  of  equations 
which  are  efficiently  solved  by  the  classical  ADI  procedures  given  by  Douglas  and 
Gunn  (Ref.  27).  The  ADI  splitting  is  usually  taken  along  coordinate  directions. 
When  this  is  the  case  mixed  derivative  terms  are  often  eliminated  on  the  implicit 
side  by  use  of  the  function  e ^  of  equation  (31). 


It  is  often  instructive  to  specialize  from  the  general  to  the  specific.  In 
so  doing  one  can  consider  the  application  of  the  scheme  given  in  Eq.  (31)  to  a 


fc'Pkx  ,  Bcpk 

system  with  =  Hi('Pi,  — £)  and  =  F  (<Pj,  -j)- 


..  This  system  covers  all  equa- 

by*'  *  d  5y* 

tions  that  do  not  have  time-like  derivatives  outside  of  the  solution  vector  itself. 
In  gas  dynamics  this  type  of  system  is  obtained  when  a  direct  solution  to  the 
steady-state  inviscid  supersonic  equations  is  desired  in  Cartesian  coordinates. 

h  u  Hj[ 


From  Eq.  (25)  one  obtains  a^  =  b^  =0,  c^ 


’  and  diJ/t  2  ’ 


and  from 


.fiH, 


Eq.  (17)  one  has  s^  =  BjiFi  =  (^ji  Fi 
one  gets 


txp. 


i' 


F. .  By  substitution  into  Eq.  (31) 


+  h  fc2Hi 


h  t>Fi 

F  -  -  - 

2  bHr  r  2 


n 


,  n+1  ,  n 

h  ■  hFi 


(34) 


with  a  choice  of  ei/t  =0.  If  one  had  considered  the  time -dependent  equations  in 
primitive  variables,  then  one  would  obtain  and  hence  the  system 


fii  . 

CFil 

£ 

/O 

IfM 

in 


^+1  -  hFin 


(35) 
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The  further  restriction  to  a  single  ordinary  differential  equation  would  lead  to 


(1  -  |  £>“  (2^^)  =  ^  «6) 


or 


„n+l  „n 


<pn+1  -  cpn  ,  „  n  5F  u  <P“~  -  qT  .. 

- - -  4  C  Fn  +  [  Fn  +  (-)  ( - — )  h]) 

h  2  C*P  h 


(37) 


n+1 


But  now  within  the  square  parenthesis  one  has  an  approximation  to  F  accurate  to 
second  order.  If  the  system  were  linear  in  F,  then  one  would  have 


<pn+1  -  cp11 


Fn  +  Fn+1 


which  is  exactly  the  Crank -Nicolson  scheme. 


(38) 
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THE  COORDINATE  SYSTEM 


The  Construction  of  Tube -like  Coordinates 

Tube-like  coordinates  will  be  constructed  in  general  to  provide  a  natural 
setting  for  the  study  of  flows  within,  between,  or  outside  of  a  set  of  prescribed 
tubes.  The  prescribed  boundary  tubes  then  become  coordinate  surfaces,  and,  as  a 
result,  the  specification  of  fluid  dynamic  boundary  conditions  is  greatly  simplified. 
Although  the  equations  of  motion  contain  more  terms  than  for  a  cartesian  system, 
this  does  not  add  excessively  to  the  run  time  of  a  program.  In  addition,  there 
must  be  some  control  over  the  resolution  of  regions  near  bounding  tubes  so  that  the 
effects  of  wall  curvature  and  the  growth  of  attached  boundary  layers  can  be  ade¬ 
quately  treated.  Such  controls  are  obtained  from  the  specification  of  coordinate 
distribution  functions  which  shall  appear  only  as  parameters  in  the  basic  geometric 
construction  of  the  coordinates.  The  basic  geometry  of  the  bounding  tubes  then 
provides  the  intrinsic  constraints  upon  the  coordinate  construction.  Since  the 
primary  goal  is  the  computation  of  fluid  flows  within  nontrivial  geometries  and  not 
the  development  of  coordinate  systems  per  se,  the  coordinates  will  be  kept  as  simple 
as  possible,  given  the  desired  generality. 

Considering  various  past  successes  of  two-dimensional  conformal  mappings  to 
obtain  coordinates,  one  might  naturally  wish  to  obtain  similar  transformations  for 
three-dimensions.  Unfortunately,  there  is  no  three-dimensional  theory  of  conformal 
transformations  analogous  to  complex  variables,  and  consequently,  in  three  dimen¬ 
sions  one  is  left  with  a  complicated  system  of  partial  differential  equations  which 
generally  would  require  numerical  solution.  To  circumvent  the  considerable  computa¬ 
tional  labor  required  for  solution  of  such  equations,  a  constructive  process  is 
used  for  the  development  of  tube-like  coordinates. 

The  first  step  in  the  construction  of  tube-like  coordinates  is  to  create  a 
suitable  family  of  two-dimensional  surfaces  which,  in  some  sense,  are  transverse  to 
a  given  centerline.  If  orthogonal  coordinates  are  desired,  then  these  surfaces  would 
have  to  bend  and  flex  as  the  tube  would  undergo  changes  in  cross  section  at  differ¬ 
ent  centerline  positions.  In  addition  to  the  problem  of  constructing  transverse 
surfaces  which  bend  and  flex,  there  is  also  the  problem  of  constructing  an  orthogonal 
grid  on  a  surface  which  has  variations  in  Gaussian  curvature,  and  hence,  is  not  flat. 
This  second  problem,  in  fact,  requires  a  more  complicated  construction  than  the 
first  which  in  itself  is  not  easy.  Thus,  the  sheer  magnitude  of  the  work  involved 
in  the  construction  of  orthogonal  coordinates  certainly  would  remove  the  desire 
for  their  use  in  fluid  dynamic  problems  which  undoubtedly  would  require  less  com¬ 
putation  in  nonorthogonal  coordinates  than  in  the  construction  of  an  orthogonal 
system  alone.  By  contrast,  if  the  transverse  surfaces  are  selected  to  be  two-dimen¬ 
sional  planes,  then  the  construction  of  coordinates  is  greatly  simplified  while  the 
fluid  dynamic  computation  is  only  marginally  different  due  to  coordinate  nonortho¬ 
gonality.  Consequently  the  coordinate  system  that  we  shall  construct  will  have 
planar  transverse  surfaces. 
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Since  each  planar  transverse  surface  is  a  linear  subspace  of  the  real 
three-dimensional  Euclidian  vector  space  R^,  any  such  plane  can  be  completely 
specified  by  any  two  spanning  linearly  independent  vectors  in  R  .  The  specifica¬ 
tion  of  the  planar  family  of  transverse  surfaces  is  then  a  result  of  a  construction 
of  two  vector  fields  along  a  given  centerline  curve  in  R  .  The  origin  of  each  plane 
is  chosen  to  coincide  with  the  associated  centerline  point.  (See  Fig.  4.)  To 
assure  that  the  planes  are  always  indeed  transverse,  it  will  be  assumed  that  they 
are  orthogonal  to  the  centerline  at  their  origins.  Ultimately,  tube-like  surfaces 
will  be  generated  by  loops  about  the  planar  origins  which  deform  in  some  way  as  we 
move  along  the  centerline  curve;  in  general,  these  tube-like  surfaces  will  not 
intersect  the  transverse  planes  orthogonally.  Thus,  only  the  centerline  direction 
determines  the  transverse  nature  of  the  cross  sectional  planes.  Specifically, 
the  centerline  tangent  vectors  form  a  vector  field  which,  at  each  point,  is  ortho¬ 
gonal  to  the  plane  of  the  two  transverse  vectors,  and  thus  each  centerline  point 
carries  a  triple  of  linearly  independent  vectors.  By  the  Gram-Schmidt  orthogonali- 
zation  procedure,  each  such  triple  of  vectors  can  be  made  into  an  orthogonal  set, 
and  hence,  an  orthonormal  set  which  is  simply  called  a  frame.  Thus,  tube-like 
coordinate  systems  are  constructed  from  a  specified  centerline  curve  and  an  asso¬ 
ciated  frame  field.  Now  the  basic  question  is  whether  there  is  a  canonical  con¬ 
struction  of  tube-like  coordinate  systems  from  either  a  given  centerline  or  a  given 
frame  field.  From  the  theory  of  space  curves  (Ref.  26),  it  is  well  known  that  for 
positive  curvature  and  specified  torsion  there  is  a  local  one-to-one  correspondence 
between  frame  fields  and  space  curves  which  pass  through  a  given  point.  Thus,  for 
nonzero  curvature,  the  centerline  space  curve  has  a  canonical  frame  field  which  is 
known  as  the  Frenet  frame.  Consequently,  the  coordinates  will  be  derived  from  the 
Frenet  frame  when  it  exists.  At  centerline  points  of  zero  curvature,  the  Frenet 
frame  is  degenerate  and  must  be  treated  specially. 

Once  the  Frenet  frame  of  the  space  curve  ~y  has  been  established,  the  unit 
normal  and  binormal  vectors  Vp  and  at  each  point  of  y  determine  a  transverse 
plane  orthogonal  to  the  unit  tangent  vector  V^.  (See  Fig.  5-)  Relative  to  any  such 
transverse  plane,  these  vectors  are  also  the  standard  orthonormal  basis.  Conse¬ 
quently,  we  can  examine  the  plane  separately  from  the  curve,  %  which  will  only 
appear  as  the  point  at  the  origin.  In  two  dimensional  functional  terminology,  the 
unit  normal  direction  can  be  considered  as  the  abscissa  and  the  unit  binormal  as  the 
ordinate;  or  more  simply,  as  x  and  y  axes,  respectively.  Since  the  tube-like 
coordinates  are  to  be  generated  from  some  family  of  tubes  encasing  the  space  curve, 
"y,  a  cross-sectional  cut  by  a  transverse  plane  produces  within  the  plane  a  family 
of  loops  about  the  origin.  We  shall  assume  that  each  loop  is  representable  by  a 
strictly  monotone  radial  function  of  angle.  In  this  regard,  a  polar  type  of  descrip 
tion  is  the  most  suitable.  But,  of  course,  the  loops  are  usually  more  complicated 
than  circles,  and  thus,  we  must  replace  the  radius  by  a  function  L  of  both  radial 
and  angular  variables  r  and  0.  Furthermore,  when  noncircular  loops  bound  a  cross 
section  of  fluid,  there  are  regions  of  varying  wall  curvature.  In  a  numerical 
solution,  it  is  desirable  to  put  proportionately  more  mesh  points  in  regions  of 
higher  curvature  than  in  regions  of  less  curvature .  Consequently ,  an  angular 
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distribution  function,  0,  is  a  good  replacement  for  the  simple  angular  specification, 

9,  of  simple  polar  coordinates.  The  net  result  is  a  generalization  where  polar 
coordinates  are  replaced  by  a  pseudo-radius,  r  and  pseudo-angle,  0.  Since  the  loops 
generally  vary  from  transverse  plane  to  transverse  plane,  the  pseudo-radii  and  angles 
must  also  be  functions  of  axial  location,  t,  on  the  centerline  space  curve,  ^(t). 

Since  the  normal  and  binormal  directions  are  usually  functions  only  of  the  centerline 
curve,  y(t),  our  loops  may  have  symmetries  that  do  not  reflect  about  either  of  these 
Frenet  directions.  Since  the  use  of  known  symmetries  is  a  great  simplification  in 
most  problems,  we  need  an  option  which  allows  one  to  define  axes  that  can  be  aligned 
in  an  optimal  way.  This  option  is  easily  established  from  the  specification  of  a 
function,  Cl(t),  which  is  a  rigid  rotation  relative  to  the  normal -binormal  directions. 

To  bring  this  development  of  tube-like  coordinates  within  the  framework  of  the 
preceding  tensor  derivations,  we  shall  use  the  notation,  y=0>  y  =r  and  y  =t  for 
pseudo-angular,  pseudo-radial,  and  axial  variables.  In  this  notation,  we  have  thus 
far  developed  (1)  a  length  factor,  L^y1,  y2,  y3),  which  is  a  generalization  of 
radius,  (2)  an  angular  distribution  function,  0=®(y1,  y3)  which  is  a  generalization 
of  angle,  (3)  a  rotation  function,  0=Q( y3),  and  (4)  the  Frenet  frame,  (V15  V2,  V^)= 

(V  (y3 ) ,  f2(y3),  Vo(y3))  upon  which  the  coordinates  are  built.  That  the  length 
factor,  L,  and  the  angular  distribution  function,  ©,  give  us  a  generalization  of 
polar  coordinates  is  obvious  since  polar  coordinates  are  easily  retrieved  by  taking 
L(yL,  y2,  y3)  =  y2  and  ©(y1,  y3)  =  y1.  It  is  also  worth  noting  that  the  angular 
distribution  function,  0,  was  chosen  to  be  independent  of  pseudo-radius,  y  . 

Although  it  is  not  immediately  evident,  we  have  removed  a  considerable  amount  of 
potential  computational  complexity  in  the  process  of  obtaining  metric  information 
by  limiting  the  number  of  derivatives  which  must  be  computed.  Furthermore,  there 
is  no  real  loss  of  flexibility  in  the  construction  of  angular  distribution  functions. 
Since  most  commonly  used  analytic  descriptions  of  loops  are,  in  fact,  controlled  by 
a  collection  of  parameters  which  depend  only  on  axial  location,  y  ,  a  knowledge  of 
only  these  parameters  is  often  sufficient  for  the  construction  of  the  angular  dis¬ 
tribution  function.  For  example,  if  the  loops  were  to  consist  of  a  family  of  con¬ 
centric  homogeneous  ellipses,  then  the  major  and  minor  axes  of  the  outermost  ellipse 
would  form  a  collection  of  two  such  parameters. 

With  the  above  functions  and  the  Frenet  frame,  the  class  of  tube-like  coordinates 
comes  directly  out  of  the  transformation 

"x  *=  "y  +  L{V2  cos  cp  +  sin  cp}  ^ 

which  transforms  curvilinear  coordinates,  y’  =  (y1,  y2 ,  y3)  into  cartesian  coordinates 
x”  =  (x1,  X2,  x3)  Where 


cpCy1,  y3)  «  ©(y1*  y^)  +  n(y3). 


(4o) 


. 
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At  each  transverse  location,  y3,  the  space  curve  vector,  y,  translates  the  origin 
to  the  space  curve.  At  a  given  pseudo-angle,  jA,  a  unit  vector,  V  cos  cp  +  V3  sin  cp,  is 
determined  by  the  sum,  cp  =  ©  +  0  of  the  radial  distribution  function,  ®,  and  the 
transverse  rotation,  Q.  This  unit  vector  sweeps  out  a  full  360  degs  in  the  trans¬ 
verse  plane  as  y1  passes  through  all  of  its  values.  Hence,  we  could  call  this  a 
direction  pointer  for  the  transverse  plane.  When  this  direction  pointer  is  scaled 
by  the  length  factor,  L,  we  obtain  a  point  of  our  transformation.  Since  the  length 
factor  depends  on  all  three  variables,  any  set  of  tube-like  surfaces  can  be  obtained 
provided,  of  course,  that  loops  are  representable  by  a  strictly  monotone  radial 
function  of  angle  and  also  that  no  two  transverse  cross  sections  are  allowed  to 

intersect. 

In  a  geometric  setting,  the  transformation  is  really  an  embedding  of  tube-like 
coordinate  systems  into  three  dimensional  Euclidian  space.  An  illustration  is 
provided  in  Fig.  5.  From  the  transformation,  it  is  also  easy  to  see  that  the  surfaces 
of  constant  y3  are  the  transverse  planes,  the  surfaces  of  constant  pseudo-angle,  y  , 
are  ruled  surfaces  generated  from  the  centerline  curve,  y ,  and  the  surfaces  of^ con¬ 
stant  pseudo-radius,  y2,  are  just  the  concentric  tubes  about  the  space  curve,  y 
Separate  illustrations  of  these  various  coordinate  surfaces  are  given  in  Figs,  6a, 

6b,  and  6c,  respectively. 


The  Length  Factor 

the  structure  of  tuhe-like  coordinates ,  the  length  factor  contains  the  in- 
Within  the  struc  specification  of  basic  geometry  and  for  the  distributional 

formation  needed  both  for  P  ^  oontrol  can  be  easily  implemented  by 

"  ff  pseudo-radial  and  angular  distr^t  ions  * 

atruction  of  the  length  factor.  In  his  w,  tM ba »  6e-e  ry  especlally 

rately  from  the  distributional  aspects.  down  into  a  number  . 

“si^e^'hounding  tubes  can  be  smoothly  generated  from  bounding  loops 
5tT~h  “it  is  sufficient  to  temporarily  =t ^alysi^ 

to  a  consideration  of  regions  be  ween tJnassumed  that  the  nondegenerate  bounding 
our  tube-like  coordinates,  we  have  impxici  y  radial  lines  emanating 

loops  do  not  pass  through  the  origin  and  are  contractable  along  radia:  li 

from  the  origin.  Thus,  we  do  not  allow  ^^^..ed  in  terms  of  polar 

more  than  once.  Consequently,  eac  oun  g  1  function  Ff  (9 )  and  the  unit 

coordinates  (r,«)  the  s^olar" description  the  contraction  process  effectively 

vector  (cos  9,  sin  6 ) .  With  P  is  a  proCess  between  the  coef- 

reduces  to  a  one-dimensional  process.  Specific  y 

ficients  of  the  unit  vector  (cos  e,  sin  .).  ‘  £ c °  j  UlV  satisfactory. 

sufficiently  the  fl»  region  is 

If  we  assume  that  no  two  tubes  join  w  nff  tubes  The  interpolation 

divisible  into  subregions  with  no  more  than  wo  *  f  6  ^  is  usuany  needed; 

process  for  two  bounding  loops  Yx,  and  y2  is,  therefore,  an 
and  is  given  by  the  simple  homotopy 


*SRS3»Ka^iP>*^^^ 
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waaswjMna* 


H  (0,r)  =  L(e,r)  (cos  e,  sin  e)  (4la) 

with  length  factor 

L(0,r)  =  rF^G)  +  (1-r)  Fg  (e)  (Mb) 

which  takes  y2(0)  =  H(0,o)  uniformly  and  smoothly  into  \1(e)  =  H(0,1)  as  r  goes  from 
0  to  1.  (See  Ref.  28.).  This  is  illustrated  in  Fig.  7.  If  the  loop  y2  is  degenerate 
then  the  coefficient  F  (0)  vanishes  and  the  length  factor  reduces  to  L(e,r)  =  rF  (0), 
we  thus  have  the  cross  section  of  a  duct  generated  by  one  loop.  By  the  continuity 
of  L  as  a  function  of  y2»  the  duct  generated  by  one  loop  y1  can  be  considered  as  a 
limit  of  annular  type  regions  between  loops  Yp  Y^_  and  y2  closes  tightly  upon  the 

origin.  This  concept  is  often  quite  useful  since  the  origin  in  coordinates  generated 
from  one  loop  suffer  from  the  same  singularity  problem  that  occurs  with  simple  polar 
coordinates.  This  singularity  can  be  circumvented,  however,  by  using  an  auxiliary 
loop  y2  which  is  near  enough  to  the  origin  to  create  a  good  approximation  to  the 
original  region.  To  preserve  overall  accuracy  in  a  numerical  computation,  ||y2II  = 
ra“cF2(0)  must  be  less  than  the  numerical  truncation  error.  In  fact,  the  well-defined 
limiting  process  would  lead  one  to  believe  that  there  would  be  no  problem  at  all  in 
taking  I  Iy2I  I  arbitrarily  small.  But  if  |  |y2I  I  is  taken  within  the  region  of  machine 
roundoff  error,  then  the  singularity  problem  may  reappear  by  default.  Consequently, 
it  is  best  to  choose  |  |  y2  1 1  to  be  much  less  than  truncation  errors  but  greater  than 
roundoff  errors. 


The  final  stage  of  length  factor  construction  is  accomplished  by  a  replacement  of 
the  polar  coordinates  r  and  6  by  radial  and  angular  distribution  functions  R(r,t)  and 
0(0  ,t)  for  axial  location  t.  Now  since  R  and  0  are  to  be  the  actual  polar  locations 
of  a  loop  we  must  reinterpret  r  and  0  as  pseudo-radial  and  pseudo-angular  locations  on 
the  same  loop.  Within  this  context  the  two-tube  length  factor  becomes 


L(6,r,t)  =  R(r,t)  Fi(0(0,t))  +  fl-R(r,t)]  F2(0(0>t))  (M) 

and  the  associated  unit  vector  becomes 

(cos  (@(8,t)  +n(t)),  sin  (0(e,t)  +  n(t))  (43) 


where  the  rotation  fl(t)  of  the  Frenet  frame  has  been  included  for  completeness. 
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The  Construction  of  Bounding  Tubes 

As  the  reader  has  just  seen,  the  construction  of  tube-like  coordinate  systems 
relies  upon  the  existence  of  smooth  bounding  tubes.  When  such  tubes  exist  at  the 
outset  of  a  problem,  the  generation  of  coordinates  is  a  straightforward  application 
of  the  development  above.  However,  if  the  bounding  tubes  are  unknown  at  the  outset, 
then  they  must  be  constructed  in  a  smooth  enough  fashion.  In  such  cases  one  is 
often  given  a  discrete  specification  of  a  sequence  of  bounding  loops  which  first 
must  be  fit  with  a  smooth  curve  and  then  must  be  joined  to  form  a  smooth  surface. 

This  circumstance  can  often  arise  out  of  the  convenience  associated  with  the  discrete 
sp.oification  of  a  surface  by  means  of  successive  cross  sectional  cuts.  If  this 
data  were  to  be  given  all  in  advance  of  the  intended  use  of  the  coordinate  system 
as  a  whole,  then  the  smoothed  cross  sectional  loops  could  be  effectively  joined  by 
fitting  them  together  witn  splines  which  can  have  interior  knots  corresponding  to 
interior  loops.  However ,  it  is  often  the  case  that  discrete  loop  data  is  only 
generated  one  station  in  advance  of  the  use  of  the  coordinate  system.  This  occurs, 
for  fiXanple,  whim  the  problism  ia  to  Jicl  v«.-  f.r  the  vlvutifl  fl  :w  fiel  *  out  tile  of  an 
ogival  body  when  the  flow  is  predominately  supersonic.  While  the  ogival  body  sur¬ 
face  is  known  in  advance,  the  location  of  the  bow  shock  is  not.  Thus  one  considers 
the  ogival  body  surface  as  the  unknown  outer  tube  which  one  wishes  to  use  for  the 
generation  of  tube-like  coordinates  to  allow  for  the  efficient  computation  of  the 
fully  viscous  flow  field.  Since  the  fl^w  field  in  a  neighborhood  of  the  bow  shock 
is  largely  inviscid,  an  inviscid  explicit  solution  is  performed  iteratively  to 
obtain  discretely  the  geometry  of  the  bow  shock  at  one  station  in  advance  of  our 
known  soltuion  and  coordinate  system.  One  is  now  left  with  a  fully  developed  bow 
shock  surface  preceeded  by  a  discrete  cross  sectional  loop  of  bow  shock  data.  The 
problem  is  to  smoothly  fit  the  loop  and  then  smoothly  join  the  result  to  obtain  a 
smooth  extension  of  the  surface.  Since  fluctuations  may  arise  from  the  discrete 
generation  of  the  tow  shook  Ista,  a  I* art,  squares  splint-  procedure  ia  to  fit 
the  loop  with  smoothness  up  to  three  continuous  derivatives.  This  type  of  least- 
squares  procedure  has  the  distinct  advantage  of  accurately  representing  the  surface 
normal  curvature  along  the  loop.  Now  one  has  a  loop  y  (0)  at  level  n+1  and  a 
bounding  tube  ending  on  a  loop  Vn(0)  at  level  n  where  SSere  are  known  derivatives 
in  the  axial  direction  t, 


One  then  attaches  a  surface  which  extends  the  tube 

The  extension 

is  accomplished  with  the  tensor  product  form 

r+1 

p  (e,t)  =  l  f.  (t)  rn+1_3  (e)  (44) 


from  7n  to  with  the  smoothness  of  three  continuous  derivatives 


which  takes  information  back  to  loop  ^  .  At  the  beginning  r  must  be  1  since  there 

are  only  two  available  loops.  The  process  continues  with  r  increasing  until  a 
desired  maximum  r  value  is  obtained.  From  there  on  r  is  assumed  to  be  constant. 
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Distribution  Functions 

When  partial  differential  equations  are  discretized  in  terms  of  differences, 
the  derivatives  are  replaced  in  some  fashion  by  difference  quotients.  A  simplifi¬ 
cation  then  leads  to  the  difference  equations  that  we  solve.  Implicitly  in  the 
discretization,  however,  is  the  assumption  that  derivatives  are  accurately  estimated 
by  secant  lines.  But  then  the  exact  solution  may  experience  drastic  variations  in 
a  short  distance.  Such  solutions  are  said  to  have  large  gradients.  In  regions 
where  the  gradients  are  large,  the  approximation  of  derivatives  by  secants  may  be 
very  poor  unless  the  particular  region  is  disected  into  smaller  regions  which  have 
reasonable  secant  approximations,  a  practice  commonly  known  as  mesh  refinement.  In 
fluid  mechanics,  the  boundary  layer  of  a  viscous  flow  around  or  through  an  object  is 
such  a  region. 

Obviously,  the  necessary  resoltuion  could  be  accomplished  by  merely  increasing 
the  number  of  points  in  a  uniform  distribution;  however,  this  would  require  excessive 
computer  time  and  storage.  Another  alternative,  known  as  the  interface  method,  is 
to  use  a  refined  mesh  only  in  the  given  region  and  then  join  it  with  the  global 
mesh.  An  improved  technique  is  to  use  coordinate  distribution  functions  which 
smoothly  distribute  mesh  points  so  that  in  some  sense  they  are  spaced  in  roughly 
an  inverse  proportion  to  the  size  of  the  gradients.  Thus,  regions  of  high  gradients 
have  proportionately  more  points  than  regions  with  smaller  gradients.  Unlike  the 
interface  method,  the  transition  between  different  mesh  lengths  is  made  continuously, 
and  as  gradually  as  possible.  Distributions  are  often  used  when  the  distributional 
transformation  is  applied  to  an  independent  variable  of  an  existing  transformation. 

The  result  is  a  new  transformation  obtained  by  composition.  With  this  approach,  the 
problem  of  mesh  point  distribution  is  replaced  by  the  problem  of  selecting  a  suitable 
set  of  distribution  functions  within  a  transformation  of  coordinates.  The  problem 
is  a  nontrivial  one  since  the  distribution  functions  should  depend  upon  the  nature 
of  the  solution  being  computed  but  are  determined  in  advance  of  the  computation. 

Thus,  some  prior  knowledge  of  the  solution  is  required.  In  flows  with  large  boundary 
layer  separation  or  with  adjacent  dissimilar  components,  the  critical  region  to  be 
resolved  is  somewhere  in  the  middle  of  the  flow.  But  the  location  of  such  regions 
is  often  unknown  at  the  outset  of  the  problem.  One  method  to  overcome  this  difficulty 
in  marching  procedures  is  to  create  the  distribution  function  at  the  next  level  based 
upon  a  knowledge  of  the  solution  at  the  present  level.  Care  must  be  taken,  however, 
to  create  a  distribution  function  that  is  sufficiently  smooth  in  the  marching 
direction. 

In  many  problems  of  practical  interest,  however,  the  regions  that  need  resolu¬ 
tion  are  known  in  advance.  Typical  examples  are  attached  boundary  layers  and  boundary 
layers  that  may  have  small  separations  or  separation  bubbles. 

Within  the  framework  of  tube-like  coordinate  systems  boundary  layer  resolution 
on  the  inner  surface  is  accomplished  by  setting 
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R(r,t)  =1  +  1  (r-l)  tanh  [D  (  ^)]  (45) 

a  x-c* 

where  d  =  d(t)  is  the  estimated  boundary  layer  growth,  or  =  a(t)  is  the  desired 
proportion  of  mesh  points  in  the  boundary  layer,  and  D  is  the  hyperbolic  damping 
factor.  The  boundary  layer  growth  d  gives  the  fraction  of  the  flow  region  occupied 
by  the  boundary  layer,  a  is  usually  taken  as  a  constant,  and  D  can  be  given  a  value 
of  about  2.  When  r  is  small,  the  radial  distribution  of  equation  (45)  reduces 
essentially  to  the  line 


(46) 


which  would  have  been  chosen  had  we  used  the  Interface  method.  As  r  approaches 
unity  the  distribution  Eq.  (45)  smoothly  approaches  unity  as  illustrated  in  Fig.  8. 


Metric  Data  for  Tube-Like  Coordinates 


The  efficient  generation  of  metric  data  is  an  important  part  of  any  solution 
procedure  involving  general  curvilinear  coordinates.  Before  a  solution  can  be 
undertaken,  the  physical  problem  must  be  specified.  Problem  specification,  however, 
involves  the  creation  of  boundary  and  initial  data  and  the  generation  of  the  equa¬ 
tions  of  motion  with  the  associated  boundary  conditions.  In  addition,  the  solution 
may  be  monitored,  examined,  or  put  under  physical  constraints.  In  all  of  these 
tasks,  the  metric  data  is  needed.  A  knowledge  of  the  metric  data  is  enough  to  com¬ 
pletely  specify  the  equations  of  motion  and  analyze  the  coordinate  invariant  direc¬ 
tions  for  the  specification  of  boundary  and  initial  conditions.  For  very  compli¬ 
cated  geometries  the  equations  of  motion  may  contain  an  inordinate  number  of  terms. 
However,  if  the  equations  are  taken  in  tensor  form,  then  the  coefficients  to  terms 
can  be  constructed  from  the  metric  data  with  the  construction  process  being  performed 
on  a  computer.  Once  a  non-trivial  term  is  constructed,  its  contribution  to  the 
desired  difference  equations  is  computed  before  searching  for  the  next  non-trivial 
term.  Sequentially,  the  process  continues  until  all  terms  in  the  equations  have 
given  their  contributions  to  the  system  of  difference  equations.  Then,  in  the  same 
fashion,  we  cycle  through  terms  in  the  boundary  conditions,  sequentially  adding  in 
their  respective  contributions.  The  result  is  the  desired  set  of  difference  equa¬ 
tions,  and  the  problem  is  effectively  reduced  to  linear  algebra.  Note  that  with 
such  methods  there  is  no  real  need  to  write  out  the  differential  equations  or  compli¬ 
cated  boundary  conditions  in  detail.  Thus,  all  we  need  to  do  is  to  generate  the 
metric  data  and  use  it. 

The  metric  data  for  tube-like  coordinates  can  be  obtained  from  the  transformation 
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sending  tube-like  coordinates  y  (y^>y^ >y3)  =  (9>r»^)  into  cartesian  coordinates 
x  =  (x-4,x2,x3)  where  cp  =  ©  +  0,  and 

y" _  y" (y3)  =  space  curve  center-line 
V2  =  Vo  (y3)  =  unit  normal  vector 

Vg  =  V3  (y^)  =  unit  binormal  vector  (48) 

L  »  L  (y1,  y2,  y3)  =  length  factor 

@  =  ©  (y1,  y3)  =  angular  distribution 

0  =  n  (y3)  *  rotation  with  respect  to  Frenet  frame. 


By  differentiation  of  the  coordinate  transformation,  we  obtain  the  Jacobian  trans¬ 
formation  which  leads  directly  to  the  transformation  rules  for  tensor  fields.  These 
rules  allow  one  to  input,  monitor,  or  extract  basic  information  from  a  solution 
procedure  involving  transformed  variables.  The  Jacobian  Transformation  is  essentially 
obtained  from  the  chain  rule  which  yields 


bx  _  &xJ  bx_  _  £x£  a 

by1  by1  bx^  by1  ^ 


(49) 


where  ui  is  the  standard  orthonormal  basis  of  constant  vector  fields,  and  ej  is  the 
naturaibasis  of  tangent  vectors  to  coordinate  curves.  With  a  slight  abuse  of 
notation,  we  have  used  "x  as  a  position  vector  in  the  definitions  of  e*  and  u^.  However, 
nothing  is  lost  since  the  covariant  derivative  of  x  ®  xJu ..  is  just  the  partial  deri¬ 
vative  of  the  x*J  summed  on  <ij.  In  terms  of  the  notation  J 


(50) 


we  have 


and  hence  the  Jacobian  matrix 


(51) 
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(®i»  ®2>  *3) 


/  9'x1 

SF 

ax2 

iF 

w 


ax1 

aF 

9x2 

SF 

9x3 

9y3 


9x1  \ 

9F  \ 

9X2 


9x3 

aF 


(52) 


Thus,  to  obtain  the  Jacobian  transformation,  we  must  calculate  the  natural  coordinate- 
wise  basis  vectors  ej. 

For  notational  convenience,  let  derivatives  of  functions  in  only  the  axial  variable 


y3  be  denoted  by  a  dot;  and  in  addition,  let 


j 


and 


ii+j+k  L 

(52a) 

Lijk  * 

(aF)* (9y2) J(9y3'K 

1 

®ij  - 

ai+d  0 

(9yl)l(9y3)J 

(52b) 

I 

gi+d  cp 

(52c) 

'Pij  * 

(9yl)i(By3)d 

1 

where  0  <  i,d,k  <  3-  By  differentiation, we  obtain 

=  (L100  cos  <P  "  1010  sin  ^)v2  +  <L100  sin  <P  +  *%0  cos  ^^3 
®2  “  (L010  cos  ^)V2  +  (L010  sin  ^  V3 

e3  =  s  +  L00]_  (V2  cos  cp  +  V3  sin  cp) 

+  L  (V2  cos  cp  +  t3  sin  tp  -  V2  cpQ1  sin  cp  +  V3  cp01  cos  cp) 


(53a) 

(53b) 


3 

j 

1 

$ 

,-v 

!■ 


(53c) 
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Since  the  pseudo- angular  vector  €  and  the  pseudo-radial  vector  e2  are  linear 
combinations  of  the  unit  normal  vector  V2  and  the  unit  binormal  vector  V3,  they  both 
lie  entirely  within  the  transverse  plane  at  axial  location  y3.  But  unlike  e-^  and 
e  ,  the  axial  coordinate  direction  vector  S3  is  not  expressed  in  terms  of  the  basic 
Frenet  frame.  Thus,  we  cannot  easily  measure  S3  relative  to  the  transverse  plane 
unless  the  y3-derivatives  of  unit  normal  and  binormal  vectors  are  expressed  in  terms 
of  the  Frenet  frame.  The  necessary  expressions  are  given  by  the  last  two  Frenet  for¬ 
mulas  (Ref.  29)  for  the  space  curve  y;  namely, 


V2  =  -  8k  V1  +  st  V3 


(54a) 


and 


(54b) 


where  s  is  arc  length,  K  is  curvature,  and  T  is  torsion. 

By  substitution,  we  obtain 

S  =  BV1  +  (LQ01  cos  cp  -  LA.  sin  cp  )?2  +  (Lqo1  sin  cp  +  LA  cos  cp  )  (55) 


sT  +  tp  and  B  =»  s  (1-kL  cos  9).  In  summation,  we  have  just  obtained  a 
basis  S  =  b^,  from  the  Frenet  frame  T^,  ?2,  ?3  to  the  natural  basis  of 

_ n  j _ j-  -  ...  ^  TJr'QTi-l  nncTir  hoctJAVOT*  'PminH  fl  f'hfl.TID’P 


where  A 

change  of  uclbao  =  u^v.  ±  ^  j 

tube-like  coordinate  vectors  e^,  t2,  £3.  Previously,  however,  we  found  a  change  of 

u_  to  the  Frenet  frame. 

3  m  * 

im  from  cartesian 


basis  ^T-  3  a®  from  the  standard  cartesian  basis  u^j  u2>  u3 
But  then  by  composition  we  obtain  the  change  of  basis  e^  =b|  aj  u, 
tangent  vectors  to  tube-like  tangent  vectors.  By  the  identificat 


cation  developed  above 


the  Jacobian  elements  are  just 


=  bJ  a“ 

aF  1  J 


(56) 


for  i,  m  =  1,2,3. 

The  metric  tensor  g.  .  is  obtained  from  the  differential  element  of  arc  length 
(ds)2  *  gijdy^dyJ.  But  alternatively,  we  have 


,dx 


(ds)*  =  dx.dx  dyi ^ =  rr-  |*j  dyidy'*  =  (e"i»^j)  dy^y^;  (57) 


39 
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and  hence,  by  linear  independence  g^  =  e . .  That  is,  the  metric  tensor  is  computed 
by  taking  the  dot  products  of  the  tube-like  tangent  vectors  e2,  This  compu¬ 

tation  is  most  efficiently  performed  when  the  are  expressed  in  terms  of  the  Frenet 
frame.  Since  the  Frenet  frame  is  orthonormal,  we  have 


gid  =  il-ej  =  (b“  Vm).(bJ  Vjj) 
Ih  particular,  we  obtain 


m  '-l  V  V  -  b®  b"?  6  .  =  b“  b' 


b“  b*  V  .V 
bi  Vm*vX 


■i 

1  j 


mi 


m 

i  i  ;  (58) 


2  2 

611  =  *100  +  ®10^  (59a) 

2 

®22  =  L010  (59b) 

g33  *  B2  +  LQOi  +  (LA)2  (59c) 

g12  c  L100  L010  (59d) 

g13  =  Ll00  L°01  +  l201°A 

823  =  L010  L001  (59f) 


which  by  symmetry  is  the  complete  list.  Note  that  the  sines  and  cosines  have  disap¬ 
peared  as  a  result  of  cross-cancellation  and  the  identity  sin2  cp  +  cos2  cp  ■  1.  Also 
note  that  the  last  three  components  listed  are  generally  nontrivial  off-diagonal 
entries.  When  any  of  these  are  non- zero,  there  is  an  angle  other  than  90°  between 
the  respective  coordinate  directions  of  the  indices.  For  example,  the  cosine  of  the 
angle  between  the  pseudo-angular  direction  e*-j_  and  the  pseudo-radial  direction  is 
given  by  the  expression 


,gI2 . 

J&11&22 


L100 


(60) 
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which  vanishes  only  when  Lq_qq  =  0*  But  when  Lj_qo  van^-s^ies  >  4  is  independent  of  angle 
and  hence,  the  loops  axe  circles.  If,  in  addition,  Q  is  a  uniform  distribution,  then 
the  cross  section  is  given  by  polar  coordinates.  In  the  same  fashion,  axially  inde¬ 
pendent  length  factors  would  lead  to  the  vanishing  of  823*  i*1  addition,  the 

distribution  and  the  rotation  do  not  change  along  a  straight  axis,  then  g13 
vanishes,  m  this  case,  the  only  non-orthogonality  that  can  occur  is  in  the  trans¬ 
verse  plane.  If  we  combine  a-H  of  the  constraints  above,  then  the  tube-like  coor¬ 
dinates  become  cylindrical  coordinates. 


The  determinant  of  the  metric  can  be  easily  obtained  from  elementary  operations 
on  determinants.  That  is, 


g  *  det  (gy)  *  LQ10 


1  *100 
*100  Lioo  + 

LqoI  ^lOO^D!  +  l2®10A 


*001 

^OC^OOl  +  Ij2®10A 

^001  +  +  ®2 


=  L, 


'010 


(M10 


l2©10a 


L201oA  (LA)2  +  Bc 


(61) 


=  (L  LqioQio0)2 
Alternatively,  we  also  have 

g  =  det  (gy)  =  det  (e^.e^) 

=  det  [(e^,e>2,e'3)t  (ei>e2’e3^ 
=  [det  (e^,e^  %)  f 


(62) 
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where  J  is  the  Jacobian  of  the  transformation.  Consequently,  J  =  L  L010®1nB*  For 
the  transformation  to  be  nonsingular  J  must  never  vanish,  and  thus  no  factor  in  J 
can  vanish.  When  L  vanishes  the  transformation  degenerates  to  the  center  line 
(Fig.  9a).  When  the  pseudo-radial  derivative  vanishes,  two  distinct  coordinate 

loops  coincide.  This  can  happen  by  an  intersection  or  a  point  of  tangency  (Fig.  9b). 
If  were  to  vanish,  then  the  angular  distribution  0  would  fail  to  have  the 
required  monotonicity  and  hence  a  loop  would  be  tangent  to  a  radial  line  (Fig.  9c). 
Finally,  if  B  vanishes,  then  two  distinct  transverse  planes  intersect  (Fig.  9^)- 
This  leads  to  the  restriction 


T  1  1 

L  cos  c?  /  -  (63) 

But  the  curve  L  cos  cp  =  ^  is  a  line  within  a  transverse  plane  (Fig.  9e).  The  line 
is  parallel  to  the  binormal  axis  and  passes  through  the  normal  axis  at  1  /k.  As  the 
transverse  planes  move  throughout  the  centerline  space  curve,  these  lines  generate 
a  ruled  surface.  Our  condition  for  nonsingularity  is  then  that  the  space  tubes 
must  never  come  into  contact  with  this  ruled  surface.  Thus,  space  tubes  must  either 
lie  below  or  beyond  the  surface.  If  they  were  to  lie  beyond  the  surface,  however, 
then  9  would  be  restricted  to  angles  between  -  90°  and  90°.  But  then  if  the  tubes 
were  generated  by  full  transverse  loops,  they  would  have  to  possess  a  bottom  and  a 
top.  This  would  imply  that  a  given  angle  would  yield  a  radial  line  with  tangents  to 
the  loop  and/or  a  multiple  intersection  with  the  loop.  (Fig.  9f.)  A  singularity 
of  this  type  is  the  same  as  a  singularity  arising  from  a  nonmonotonic  angular  distri¬ 
bution.  Consequently,  the  restriction  is  that  the  space  tubes  must  lie  below  the 
ruled  surface,  or  algegraically  that 


L  cos  <p 


1 

K 


(64) 


The  space-tubes  are  then  assumed  to  encase  the  space  curve  y  and  never  come  into 
contact  with  the  above  ruled  surface.  For  tubes  that  have  convex  loops  symmetrically 
centered  about  the  space  curve  centerline,  the  restriction  reduces  to  the  statement 
that  the  length  factor  must  be  bounded  in  the  normal  direction  by  the  radius  of  curva¬ 
ture  (Fig.  9g).  In  general,  if  a  tube  fails  to  meet  the  inequallity  condition 
(Eq.  64),  then  it  can  be  pushed  off  center  until  this  condition  is  met.  However,  if 
the  condition  is  so  severe  that  the  centerline  must  be  taken  extremely  close  to  the 
normal  side  of  a  bounding  tube,  then  coordinates  will  become  unnaturally  spaced. 

Since  each  loop  (about  this  offset  origin)  intersects  the  normal  line  once  on  each 
side  of  the  origin,  the  mesh  spacing  drastically  changes.  This  would  tend  to  cause 
an  over-resolution  on  the  normal  side  of  the  origin  and  an  under-resolution  on  the 
opposite  side.  In  such  a  case,  it  is  best  to  generate  the  same  bounding  tube  but 
with  a  new  centerline  with  smaller  curvature.  This  is  always  possible  since  we  must 
only  offset  the  old  centerline  in  a  convex  direction  which  expands  the  curve  rather 
than  shrinking  it.  This  is  illustrated  in  Fig.  10. 
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If  there  are  no  coordinate  singularities,  then  we  can  easily  obtain  the  inverse 
of  the  metric.  In  a  sequential  order  we  have 


s33  -  5,  (65a) 

g13  =  -  £  g33 

®L0  8 


(65d) 
(65e) 

■”  ■  1*W  *  W  W'  •  i 

The  e  •-direction  covariant  derivative  Dj  of  the  vector  ei  is  again  a  vector  and 
hence  is  esqpressible  in  terms  of  the  same  basis  e^,  eg,  e^.  Specifically,  we  have 


(65b) 

(65c) 


m 

D.e .  =  r.  .  e  (66) 

i  J  ij  m  v 

m 

Hhere  the  coefficients  ^  are  known  as  Chris  toff  el  symbols.  This  covariant  derivative 
measures  the  rate  of  change  of  e.  along  a  coordinate  curve  in  the  direction  of  e^ 

This  coordinate  curve  is  an  inteiral  curve  of  e^.  It  is  obtained  by  fixing  all  except 
the  i^h  variable  in  the  transformation.  We  shall  assume  that  the  covariant  derivative 
is  the  natural  one  derivable  from  the  metric.  This  is  known  as  the  Levi-Civita 
connection  (Ref.  29).  The  Christoff el  symbols  for  this  covariant  derivative  are 
given  by  the  formula 
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{  + 
dyi  SyJ 


deij 

ay® 


} 


(67) 


This  formula  is  easily  obtained  by  differentiating  a  ®i‘^j  with  respect  to  y®1, 
permuting  all  three  of  these  indices,  forming  the  sum  in  parenthesis,  applying  sym¬ 
metry  to  the  lower  indices  of  the  Christoffel  symbols,  and  then  applying  the  inverse 
metric.  With  some  calculation,  we  obtain  the  non- zero  Christoffel  symbols  directly 
from  the  above  formula.  For  notational  convenience,  let 


B  -  di+j+kg 

iJk  "  (ayl)iOy2)J(ay3)k 


(68) 


for  0<i,j,k  <  2.  Then  we  have 


r3  _  B100 
13  ~ 


r1  =  ^io 
21  — 


r3  =  Boio 


Loio 


i_ 

®10 


r  1  _  2  LiQQ  +  ^20 
11  =  L  01O 


r  ]  =  1001  +  ©ii  +  a__  (iooo)  -  r  3  ) 

31  1  515  ®10  1  « 


r  2  ,  1101  .  WE  r  1  +  f_  r  3  ]  .  +  g23  B  B 

13  kl0  kilo  31  Q10  13  L010 
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2  JL_  ,  1  2V 

ril  =  L010  ^00  "  L100  rii  "  »L0> 


„  2  _i ,  _  .2  L1Q0  , :  2  Lnm  Av  , 

r33  - 1010 1  bos  -  “ -  <A  +  ■  > ) 


(70) 


+  (  -  s21  hoo  -  e22  B010  *  g23  Boon.)  B 


33  ®io 


(A  +  — ^221 —  )  +  (_  g11  B100  -  g12  Bq10  +  g13  B 


"OOP 


and 


r33  *  (  “  831  B100  -  g32  BQ10  +  g33  B^) 


B 


When  viscous  calculations  are  to  be  done,  we  also  need  the  derivatives  of  the  above 
Christoff el  symbols.  This  is  a  straight  forward  but  tedious  process.  The  deriva¬ 
tive  expressions  are  rather  lengthy  and  we  shall,  therefore,  not  enumerate  thera.^  One 
may  also  note  that  the  independence  of  0  from  r  and  of  R  from  0  leads  to  a  detinite 
simplification  in  the  Christoff el  sumbols.  Had  this  independence  not  existed,  then 
the  above  calculations  would  have  been  even  more  complex. 


The  geometric  interpretation  of  the  Christoffel  symbols  above  is  quite  natural. 
For  example,  consider  the  special  case  when  the  centerline  space  curve  ^  is  a^straight 
line.  Since  the  contravariant  vector  e3  is  defined  by  the  relation  e  *ei  *  6.,  it 
is  perpendicular  to  the  transverse  plane.  Consequently,  e*3  is  a  constant  vector 
which  is  parallel  to  the  straight  centerline.  But  then  the  covariant  derivative  of 
a  constant  vector  vanishes.  By  the  Liebnitz  rule  for  covariant  derivatives,  we  then 

have 


0  .  Di  <»3)  .  D^.eS) 


=  (Di'ejJ.'e3  +  ej.(D£e3) 


=  (I'ij  +  ® 


6 


3 


m 


(71) 
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for  all  i,  j  =>  1,2,3.  This  fact  can  also  be  verified  from  our  list  of  analytic 
expressions  of  Christoffel  symbols.  The  geometric  interpretation  is  that  the  covari¬ 
ant  derivatives  of  natural  basis  vectors  ?1,  eg,  e^  have  no  axial  projection  and 
hence  are  transverse.  In  summary,  we  have  just  shown  that  with  straight  center- 
lines,  no  vector  from  the  natural  basis  ?2,  can  change  in  the  axial  direction. 
The  only  such  vector  with  a  non- trivial  axial  projection  is,  however,  the  natural 
vector  e,  to  axially  generated  coordinate  curves.  Consequently,  the  axial  projection 
of  e^  along  such  a  coordinate  curve  is  a  constant  as  illustrated  in  Fig.  11. 
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TREATMENT  OF  THE  BOW  SHOCK 


Although  the  shape  of  the  body  is  specified  at  the  start  of  the  problem,  the 
shape  of  the  bow  shock  cannot  be  specified  until  effects  of  the  body  on  the  flow 
within  the  bow  shock  are  evaluated.  The  shock  shape  is  therefore  evaluated  at  each 
new  marching  station  based  on  information  already  computed.  The  shock  surface  is 
adopted  as  the  outer  coordinate  surface  and  is  used  to  determine  the  necessary 
metric  information  for  the  tube-like  coordinates.  The  governing  equations  are  then 
solved  in  the  annular  region  between  the  ogival  body  and  the  bow  shock  by  marching 
from  one  transverse  plane  to  the  next,  proceeding  in  the  nominally  streamwise 
direction. 


The  bow  shock  is  computed  as  a  discontinuity  satisfying  the  classical  Rankine- 
Hugoniot  relations.  The  intersection  of  this  shock  and  a  transverse  computational 
plane  is  a  loop  represented  by  discrete  grid  points.  Provided  that  a  given  grid 
point  on  this  loop  is  outside  the  "zone  of  influence"  of  the  neighboring  points  on 
the  loop,  the  shock  solution  at  the  given  point  is  independent  of  the  solution  at 
adjacent  points.  This  "zone  of  influence"  assumption  is  valid  over  a  wide  range 
of  flow  conditions  and  consequently  is  not  a  limiting  assumption.  Thus  the  shock 
radius  Yn+i(6)  at  each  point  in, the  n+1  transverse  plane  can  be  evaluated  inde¬ 
pendently  by  a  pointwise  iteration  procedure. 

The  iteration  at  each  circumferential  location  in  the  n+1  plane  proceeds  by 
first  locally  extending  the  shock  surface  from  the  most  recently  evaluated  trans¬ 
verse  computational  plane,  n>  to  a  point  in  the  n+1  plane.  The  extension  of  the 
shock  surface  includes  the  point  being  evaluated  in  the  n+1  plane,  but  does  not 
extend  circumferentially  to  the  neighboring  points.  This  extension  is  a  first  guess 
for  the  shock  location  at  a  circumferential  point  and  hence  for  a  point  on  the  outer 
tube-like  coordinate  surface  given  by  Yn+i(e)  in  Eq.  44. 


Given  a  guess  at  the  shock  location,  the  axial  mass  flux  inside  the  shock  can 
be  computed  by  two  methods.  First,  an  application  of  the  Rankine-Hugoniot  condi¬ 
tions  produces  a  value  of  the  axial  mass  flux  based  only  on  the  shock  shape  and  the 
flow  properties  outside  the  shock.  Second,  an  application  of  a  compatibility  condi¬ 
tion  produces  a  second  value  of  this  flux  that  depends  only  on  the  shock  shape  and 
the  flow  properties  inside  the  shock.  The  shock  location  is  then  adjusted  iteratively 
until  the  axial  mass  flux  inside  the  shock  computed  by  the  two  methods  is  the  same.. 
This  iteration  for  the  shock  location  is  repeated  at  each  of  the  circumferential  grid 
points  to  produce  a  ring  of  discrete  points  at  the  n+1  station  which  collectively 
determine  the  shock  surface.  These  discrete  points  must  be  represented  by  a  con- 
tinuous  smooth  curve  to  provide  the  information  required  to  construct  the  coordinate 
system.  For  this  purpose  a  least  squares -spline  curve  fitting  routine  is  employed. 
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FIGURE  6a.  TRANSVERSE  PLANAR  CUTS  OF  CONSTANT  AXIAL  LOCATION  y3 


FIGURE  6b.  RULED  SURFACE  OF  CONSTANT  PSEUDO-ANGLE  y1 


FIGURE  6c.  TUBE-LIKE  SURFACES  OF  CONSTANT  PSEUDO-RADIUS  y2 
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FIGURE  7.  LINEAR  HOMOTOPY  BETWEEN  TWO  LOOPS 
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FIGURE  11.  EQUAL  PROJECTIONS 
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